Fast MC Tool (Ruben)
authoramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Apr 2012 08:11:28 +0000 (08:11 +0000)
committeramastros <amastros@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 11 Apr 2012 08:11:28 +0000 (08:11 +0000)
ITS/UPGRADE/KMCDetector.cxx [new file with mode: 0755]
ITS/UPGRADE/KMCDetector.h [new file with mode: 0755]
ITS/UPGRADE/PrepSummaryKMC.C [new file with mode: 0755]
ITS/UPGRADE/prepPlots1KMC.C [new file with mode: 0755]
ITS/UPGRADE/readmeKMC [new file with mode: 0644]
ITS/UPGRADE/testDetKMC.C [new file with mode: 0644]

diff --git a/ITS/UPGRADE/KMCDetector.cxx b/ITS/UPGRADE/KMCDetector.cxx
new file mode 100755 (executable)
index 0000000..8879141
--- /dev/null
@@ -0,0 +1,2916 @@
+#include "KMCDetector.h"\r
+#include <TMath.h>\r
+#include <TH1F.h>\r
+#include <TProfile.h>\r
+#include <TF1.h>\r
+#include <TMatrixD.h>\r
+#include <TGraph.h>\r
+#include <TAxis.h>\r
+#include <TFormula.h>\r
+#include <TCanvas.h>\r
+#include <TEllipse.h>\r
+#include <TText.h>\r
+#include <TRandom3.h>\r
+#include <TGraphErrors.h>\r
+#include <TStopwatch.h>\r
+\r
+/***********************************************************\r
+\r
+Fast Simulation tool for Inner Tracker Systems\r
+\r
+***********************************************************/\r
+\r
+\r
+#define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector\r
+\r
+#define Luminosity    1.e27       // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )\r
+#define SigmaD        6.0         // Size of the interaction diamond (cm) (LHC = 6.0 cm)\r
+#define dNdEtaMinB    1//950//660//950           // Multiplicity per unit Eta  (AuAu MinBias = 170, Central = 700)\r
+// #define dNdEtaCent    2300//15000 //1600//2300        // Multiplicity per unit Eta  (LHC at 5.5 TeV not known)\r
+\r
+#define CrossSectionMinB         8    // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)\r
+#define AcceptanceOfTpcAndSi     1 //1//0.60 //0.35  // Assumed geometric acceptance (efficiency) of the TPC and Si detectors\r
+#define UPCBackgroundMultiplier  1.0   // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)\r
+#define OtherBackground          0.0   // Increase multiplicity in detector (0.0 to 1.0 * minBias)  (eg 0.0)\r
+#define EfficiencySearchFlag     2     // Define search method:\r
+                                       // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.  \r
+\r
+#define PionMass                 0.139  // Mass of the Pion\r
+#define KaonMass                 0.498  // Mass of the Kaon\r
+#define D0Mass                   1.865  // Mass of the D0\r
+\r
+//TMatrixD *probKomb; // table for efficiency kombinatorics\r
+\r
+ClassImp(KMCProbe)\r
+\r
+Int_t    KMCProbe::fgNITSLayers = 0;\r
+Double_t KMCProbe::fgMissingHitPenalty = 2.;\r
+//__________________________________________________________________________\r
+KMCProbe& KMCProbe::operator=(const KMCProbe& src) \r
+{\r
+  if (this!=&src) {\r
+    AliExternalTrackParam::operator=(src);\r
+    fMass = src.fMass;\r
+    fChi2 = src.fChi2;\r
+    fHits = src.fHits;\r
+    fFakes = src.fFakes;\r
+    fNHits = src.fNHits;\r
+    fNHitsITS = src.fNHitsITS;\r
+    fNHitsITSFake = src.fNHitsITSFake;\r
+    fInnLrCheck     = src.fInnLrCheck;\r
+    for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];\r
+  }\r
+  return *this;\r
+}\r
+\r
+//__________________________________________________________________________\r
+KMCProbe::KMCProbe(KMCProbe& src) \r
+  : AliExternalTrackParam(src),\r
+    fMass(src.fMass),\r
+    fChi2(src.fChi2),\r
+    fHits(src.fHits),\r
+    fFakes(src.fFakes),\r
+    fNHits(src.fNHits),\r
+    fNHitsITS(src.fNHitsITS),\r
+    fNHitsITSFake(src.fNHitsITSFake),\r
+    fInnLrCheck(src.fInnLrCheck)\r
+{\r
+  for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];\r
+}\r
+\r
+//__________________________________________________________________________\r
+void KMCProbe::ResetCovMat()\r
+{\r
+  // reset errors\r
+  double *trCov  = (double*)GetCovariance();\r
+  double *trPars = (double*)GetParameter();\r
+  const double kLargeErr2Coord = 50*50;\r
+  const double kLargeErr2Dir = 0.6*0.6;\r
+  const double kLargeErr2PtI = 0.5*0.5;\r
+  for (int ic=15;ic--;) trCov[ic] = 0.;\r
+  trCov[kY2]   = trCov[kZ2]   = kLargeErr2Coord; \r
+  trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;\r
+  trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];\r
+  //\r
+}\r
+\r
+//__________________________________________________________________________\r
+void KMCProbe::Print(Option_t* option) const\r
+{\r
+  printf("M=%.3f Chi2=%7.2f (Norm:%6.2f|%d) Hits: Total:%d ITS:%d ITSFakes:%d | Y:%+8.4f Z: %+8.4f |", \r
+        fMass,fChi2,GetNormChi2(kTRUE),fInnLrCheck,fNHits,fNHitsITS,fNHitsITSFake, GetY(),GetZ());\r
+  for (int i=0;i<fgNITSLayers;i++) {\r
+    if (!IsHit(i)) printf(".");\r
+    else printf("%c",IsHitFake(i) ? '-':'+');\r
+  }\r
+  printf("|%s\n",IsKilled() ? " KILLED":"");\r
+  TString opt = option;\r
+  if (!opt.IsNull()) AliExternalTrackParam::Print(option);\r
+}\r
+\r
+//__________________________________________________________________________\r
+Int_t KMCProbe::Compare(const TObject* obj) const\r
+{\r
+  // compare to tracks\r
+  const KMCProbe* trc = (KMCProbe*) obj;\r
+  if (trc->IsKilled()) {\r
+    if (IsKilled()) return 0;\r
+    return -1;\r
+  }\r
+  else if (IsKilled()) return 1;\r
+  double chi2a = GetNormChi2(kTRUE);\r
+  double chi2b = trc->GetNormChi2(kTRUE);\r
+  if (chi2a<chi2b) return  -1;\r
+  if (chi2a>chi2b) return   1;\r
+  return 0;\r
+}\r
+\r
+//__________________________________________________________________________\r
+Bool_t KMCProbe::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const\r
+{\r
+  // Get local X of the track position estimated at the radius lab radius r. \r
+  // The track curvature is accounted exactly\r
+  //\r
+  // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)\r
+  // 0  - take the intersection closest to the current track position\r
+  // >0 - go along the track (increasing fX)\r
+  // <0 - go backward (decreasing fX)\r
+  //\r
+  // special case of R=0\r
+  if (r<kAlmost0) {x=0; return kTRUE;}\r
+\r
+  const double* pars = GetParameter();\r
+  const Double_t &fy=pars[0], &sn = pars[2];\r
+  //\r
+  double fx = GetX();\r
+  double crv = GetC(bz);\r
+  if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track\r
+    if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis\r
+      double det = (r-fx)*(r+fx);\r
+      if (det<0) return kFALSE;     // does not reach raduis r\r
+      x = fx;\r
+      if (dir==0) return kTRUE;\r
+      det = TMath::Sqrt(det);\r
+      if (dir>0) {                       // along the track direction\r
+       if (sn>0) {if (fy>det)  return kFALSE;} // track is along Y axis and above the circle\r
+       else      {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle\r
+      }\r
+      else if(dir>0) {                                    // agains track direction\r
+       if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis\r
+        else if (fy>det)  return kFALSE;        // track is against Y axis\r
+      }\r
+    }\r
+    else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis\r
+      double det = (r-fy)*(r+fy);\r
+      if (det<0) return kFALSE;     // does not reach raduis r\r
+      det = TMath::Sqrt(det);\r
+      if (!dir) {\r
+       x = fx>0  ? det : -det;    // choose the solution requiring the smalest step\r
+       return kTRUE;\r
+      }\r
+      else if (dir>0) {                    // along the track direction\r
+       if      (fx > det) return kFALSE;  // current point is in on the right from the circle\r
+       else if (fx <-det) x = -det;       // on the left\r
+       else               x =  det;       // within the circle\r
+      }\r
+      else {                               // against the track direction\r
+       if      (fx <-det) return kFALSE;  \r
+       else if (fx > det) x =  det;\r
+       else               x = -det;\r
+      }\r
+    }\r
+    else {                                 // general case of straight line\r
+      double cs = TMath::Sqrt((1-sn)*(1+sn));\r
+      double xsyc = fx*sn-fy*cs;\r
+      double det = (r-xsyc)*(r+xsyc);\r
+      if (det<0) return kFALSE;    // does not reach raduis r\r
+      det = TMath::Sqrt(det);\r
+      double xcys = fx*cs+fy*sn;\r
+      double t = -xcys;\r
+      if (dir==0) t += t>0 ? -det:det;  // chose the solution requiring the smalest step\r
+      else if (dir>0) {                 // go in increasing fX direction. ( t+-det > 0)\r
+       if (t>=-det) t += -det;         // take minimal step giving t>0\r
+       else return kFALSE;             // both solutions have negative t\r
+      }\r
+      else {                            // go in increasing fx direction. (t+-det < 0)\r
+       if (t<det) t -= det;            // take minimal step giving t<0\r
+       else return kFALSE;             // both solutions have positive t\r
+      }\r
+      x = fx + cs*t;\r
+    }\r
+  }\r
+  else {                                 // helix\r
+    // get center of the track circle\r
+    double tR = 1./crv;   // track radius (for the moment signed)\r
+    double cs = TMath::Sqrt((1-sn)*(1+sn));\r
+    double x0 = fx - sn*tR;\r
+    double y0 = fy + cs*tR;\r
+    double r0 = TMath::Sqrt(x0*x0+y0*y0);\r
+    //    printf("Xc:%+e Yc:%+e\n",x0,y0);\r
+    //\r
+    if (r0<=kAlmost0) {\r
+      AliDebug(2,Form("r0 = %f",r0));\r
+      return kFALSE;\r
+    }            // the track is concentric to circle\r
+    tR = TMath::Abs(tR);\r
+    double tR2r0 = tR/r0;\r
+    double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);\r
+    double det = (1.-g)*(1.+g);\r
+    if (det<0) {\r
+      AliDebug(2,Form("g=%f tR=%f r0=%f\n",g,tR, r0));\r
+      return kFALSE;\r
+    }         // does not reach raduis r\r
+    det = TMath::Sqrt(det);\r
+    //\r
+    // the intersection happens in 2 points: {x0+tR*C,y0+tR*S} \r
+    // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det\r
+    // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)\r
+    //\r
+    double tmp = 1.+g*tR2r0;\r
+    x = x0*tmp; \r
+    double y = y0*tmp;\r
+    if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique\r
+      double dfx = tR2r0*TMath::Abs(y0)*det;\r
+      double dfy = tR2r0*x0*TMath::Sign(det,y0);\r
+      if (dir==0) {                    // chose the one which corresponds to smallest step \r
+       double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0\r
+       if (delta<0) x += dfx;\r
+       else         x -= dfx;\r
+      }\r
+      else if (dir>0) {  // along track direction: x must be > fx\r
+       x -= dfx; // try the smallest step (dfx is positive)\r
+       if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;\r
+      }\r
+      else { // backward: x must be < fx\r
+       x += dfx; // try the smallest step (dfx is positive)\r
+       if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;\r
+      }\r
+    }\r
+    else { // special case: track touching the circle just in 1 point\r
+      if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE; \r
+    }\r
+  }\r
+  //\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________\r
+Bool_t KMCProbe::PropagateToR(double r, double b, int dir) \r
+{\r
+  // go to radius R\r
+  //\r
+  double xR = 0;\r
+  double rr = r*r;\r
+  int iter = 0;\r
+  const double kTiny = 1e-4;\r
+  while(1) {\r
+    if (!GetXatLabR(r ,xR, b, dir)) {\r
+      //      printf("Track with pt=%f cannot reach radius %f\n",Pt(),r);\r
+      //      Print("l");\r
+      return kFALSE;\r
+    }\r
+    \r
+    if (!PropagateTo(xR, b)) {\r
+      if (AliLog::GetGlobalDebugLevel()>2) {\r
+       printf("Failed to propagate to X=%f for R=%f\n",xR,r); \r
+       Print("l"); \r
+      }\r
+      return kFALSE;\r
+    }\r
+    double rcurr2 = xR*xR + GetY()*GetY();\r
+    if (TMath::Abs(rcurr2-rr)<kTiny || rr<kAlmost0) return kTRUE;\r
+    //\r
+    // two radii correspond to this X...\r
+    double pos[3]; GetXYZ(pos);\r
+    double phi = TMath::ATan2(pos[1],pos[0]);\r
+    if (!Rotate(phi)) {\r
+      if (AliLog::GetGlobalDebugLevel()>2) {\r
+       printf("Failed to rotate to %f to propagate to R=%f\n",phi,r); \r
+       Print("l"); \r
+      }\r
+      return kFALSE;\r
+    }\r
+    if (++iter>8) {\r
+      if (AliLog::GetGlobalDebugLevel()>2) {\r
+       printf("Failed to propagate to R=%f after %d steps\n",r,iter); \r
+       Print("l"); \r
+      }\r
+      return kFALSE;\r
+    }\r
+  } \r
+  return kTRUE;\r
+}\r
+\r
+\r
+//__________________________________________________________________________\r
+Bool_t KMCProbe::CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward)\r
+{\r
+  //  printf("before at r=%.1f p=%.4f\n",lr->fR, P());\r
+  if (AliExternalTrackParam::CorrectForMeanMaterial(lr->fx2X0, inward ? lr->fXRho : -lr->fXRho, GetMass() , kTRUE)) {\r
+    //  printf("after  at r=%.1f p=%.4f\n",lr->fR, P());\r
+    return kTRUE;\r
+  }\r
+  AliDebug(2,Form("Failed to apply material correction, X/X0=%.4f", lr->fx2X0));\r
+  if (AliLog::GetGlobalDebugLevel()>1) Print();\r
+  return kFALSE;\r
+}\r
+\r
+/////////////////////////////////////////////////////////////////////////////\r
+ClassImp(KMCCluster)\r
+\r
+//_________________________________________________________________________\r
+KMCCluster::KMCCluster(KMCCluster &src) \r
+: TObject(src),\r
+  fY(src.fY),fZ(src.fZ),fX(src.fX),fPhi(src.fPhi)\r
+{}\r
+\r
+//__________________________________________________________________________\r
+KMCCluster& KMCCluster::operator=(const KMCCluster& src) \r
+{\r
+  if (this!=&src) {\r
+    TObject::operator=(src);\r
+    fY = src.fY;\r
+    fZ = src.fZ;\r
+    fX = src.fX;\r
+    fPhi = src.fPhi;\r
+  }\r
+  return *this;\r
+}\r
+\r
+//_________________________________________________________________________\r
+void KMCCluster::Print(Option_t *) const \r
+{\r
+  printf(" Local YZ = (%3.4lf,%3.4lf) | X=%3.4lf  phi: %+.3f %s\n",fY,fZ,fX,fPhi,IsKilled()?"Killed":""); \r
+}\r
+\r
+/////////////////////////////////////////////////////////////////////////////\r
+ClassImp(KMCLayer)\r
+\r
+Double_t KMCLayer::fgDefEff = 1.0;\r
+//__________________________________________________________________________\r
+KMCLayer::KMCLayer(char *name) : \r
+  TNamed(name,name),fR(0),fx2X0(0),fPhiRes(0),fZRes(0),fEff(0),fIsDead(kFALSE),fType(-1),fActiveID(-1),fSig2EstD(999),fSig2EstZ(999),\r
+  fClCorr(),fClMC(),fClBg("KMCCluster",5), fTrCorr(), fTrMC("KMCProbe",5)\r
+{\r
+  Reset();\r
+}\r
+\r
+//__________________________________________________________________________\r
+void KMCLayer::Reset() \r
+{\r
+  fTrCorr.Reset();\r
+  fClCorr.Reset();\r
+  ResetMC();\r
+  fSig2EstD = fSig2EstZ = 999;\r
+  //\r
+}\r
+\r
+//__________________________________________________________________________\r
+KMCProbe* KMCLayer::AddMCTrack(KMCProbe* src) \r
+{\r
+  int ntr = GetNMCTracks(); \r
+  KMCProbe* prb = 0;\r
+  if (src) prb = new(fTrMC[ntr]) KMCProbe(*src);\r
+  else     prb = new(fTrMC[ntr]) KMCProbe();\r
+  if (!IsDead()) prb->ResetHit(GetActiveID());\r
+  return prb;\r
+}\r
+\r
+//__________________________________________________________________________\r
+void KMCLayer::Print(Option_t *opt) const\r
+{\r
+  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);\r
+  TString opts = opt; opts.ToLower();\r
+  if (opts.Contains("c")) {\r
+    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());\r
+  }\r
+}\r
+\r
+/////////////////////////////////////////////////////////////////////////////\r
+Double_t KMCDetector::fgVtxConstraint[2]={-1,-1};\r
+\r
+ClassImp(KMCDetector)\r
+KMCDetector::KMCDetector() :\r
+TNamed("test_detector","detector"),\r
+  fNLayers(0),\r
+  fNActiveLayers(0),\r
+  fNActiveITSLayers(0),\r
+  fLastActiveLayer(-1),\r
+  fLastActiveITSLayer(-1),\r
+  fLastActiveLayerTracked(-1),\r
+  fBFieldG(5.),\r
+  fLhcUPCscale(1.0),    \r
+  fIntegrationTime(0.02), // in ms\r
+  fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence\r
+  fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle\r
+  fParticleMass(0.140),    // Standard: pion mass \r
+  fMaxRadiusSlowDet(10.),\r
+  fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers\r
+  fAtLeastFake(1),       // if at least x fakes, track is considered fake ...\r
+  fdNdEtaCent(2300),\r
+  fMaxChi2Cl(25.),\r
+  fMaxNormChi2NDF(5.),\r
+  fMinITSHits(4),\r
+//\r
+  fMaxChi2ClSQ(4.), // precalulated internally\r
+  fMaxSeedToPropagate(300),\r
+//\r
+  fUseBackground(kFALSE),\r
+  fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),\r
+  fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),\r
+  fDensFactorEta(1.),\r
+//\r
+  fUpdCalls(0),\r
+  fHMCLrResidRPhi(0),\r
+  fHMCLrResidZ(0),\r
+  fHMCLrChi2(0),\r
+  fPattITS(0)\r
+{\r
+  //\r
+  // default constructor\r
+  //\r
+  //  fLayers = new TObjArray();\r
+  RequireMaxChi2Cl(fMaxChi2Cl); // just to precalulate default square\r
+}\r
+\r
+KMCDetector::KMCDetector(char *name, char *title)\r
+  : TNamed(name,title),\r
+    fNLayers(0),\r
+    fNActiveLayers(0),\r
+    fNActiveITSLayers(0),\r
+    fLastActiveLayer(-1),\r
+    fLastActiveITSLayer(-1),    \r
+    fLastActiveLayerTracked(-1),\r
+    fBFieldG(5.0),\r
+    fLhcUPCscale(1.0),\r
+    fIntegrationTime(0.02),  // in ms\r
+    fConfLevel(0.0027),      // 0.27 % -> 3 sigma confidence\r
+    fAvgRapidity(0.45),      // Avg rapidity, MCS calc is a function of crossing angle\r
+    fParticleMass(0.140),     // Standard: pion mass\r
+    fMaxRadiusSlowDet(10.),\r
+    fAtLeastCorr(-1),     // if -1, then correct hit on all ITS layers\r
+    fAtLeastFake(1),       // if at least x fakes, track is considered fake ...\r
+    fdNdEtaCent(2300),\r
+    fMaxChi2Cl(9.),\r
+    fMaxNormChi2NDF(5.),\r
+    fMinITSHits(4),\r
+    //\r
+    fMaxChi2ClSQ(3.), // precalulated internally\r
+    fMaxSeedToPropagate(50),\r
+    //\r
+    fUseBackground(kFALSE),\r
+    fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),\r
+    fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),\r
+    fDensFactorEta(1.),\r
+    //\r
+    fUpdCalls(0),\r
+    fHMCLrResidRPhi(0),\r
+    fHMCLrResidZ(0),\r
+    fHMCLrChi2(0),\r
+    fPattITS(0)\r
+{\r
+  //\r
+  // default constructor, that set the name and title\r
+  //\r
+  //  fLayers = new TObjArray();\r
+}\r
+KMCDetector::~KMCDetector() { // \r
+  // virtual destructor\r
+  //\r
+  //  delete fLayers;\r
+}\r
+\r
+void KMCDetector::InitMCWatchHistos()\r
+{\r
+  // init utility histos used for MC tuning\r
+  enum {kNBinRes=1000};\r
+  const double kMaxResidRPhi=1.0,kMaxResidZ=1.0,kMaxChi2=50;\r
+  int nlr = fNActiveITSLayers;\r
+  TString nam = "mc_residrhi";\r
+  TString tit = "MC $Delta Cl-Tr R#phi";\r
+  fHMCLrResidRPhi = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidRPhi,nlr,-0.5,nlr-0.5);\r
+  fHMCLrResidRPhi->GetXaxis()->SetTitle("cl-tr #Delta r#phi");\r
+  fHMCLrResidRPhi->Sumw2();\r
+    //\r
+  nam = "mc_residz";\r
+  tit = "MC $Delta Cl-Tr Z";\r
+  fHMCLrResidZ = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidZ,nlr,-0.5,nlr-0.5);\r
+  fHMCLrResidZ->GetXaxis()->SetTitle("cl-tr #Delta z");\r
+  fHMCLrResidZ->Sumw2();\r
+    //\r
+  nam = "mc_chi2";\r
+  tit = "MC #chi^{2} Cl-Tr Z";\r
+  fHMCLrChi2 = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxChi2,nlr,-0.5,nlr-0.5);\r
+  fHMCLrChi2->GetXaxis()->SetTitle("cl-tr #chi^{2}");\r
+  fHMCLrChi2->Sumw2();\r
+  //\r
+  SetBit(kUtilHisto);\r
+}\r
+\r
+\r
+void KMCDetector::AddLayer(char *name, Float_t radius, Float_t x2X0, Float_t xrho, Float_t phiRes, Float_t zRes, Float_t eff) {\r
+  //\r
+  // Add additional layer to the list of layers (ordered by radius)\r
+  // \r
+\r
+  KMCLayer *newLayer = (KMCLayer*) fLayers.FindObject(name);\r
+\r
+  if (!newLayer) {\r
+    newLayer = new KMCLayer(name);\r
+    newLayer->fR = radius;\r
+    newLayer->fx2X0 = x2X0;\r
+    newLayer->fXRho  = xrho;\r
+    newLayer->fPhiRes = phiRes;\r
+    newLayer->fZRes = zRes;\r
+    eff = TMath::Min(1.f,eff);\r
+    newLayer->fEff = eff <0 ? KMCLayer::GetDefEff() : eff;\r
+    newLayer->fActiveID = -2;\r
+    TString lname = name;\r
+    newLayer->fType = KMCLayer::kTypeNA;\r
+    if      (lname.Contains("tpc")) newLayer->fType = KMCLayer::kTPC;\r
+    else if (lname.Contains("its")) newLayer->fType = KMCLayer::kITS;\r
+    if (lname.Contains("vertex")) newLayer->SetBit(KMCLayer::kBitVertex);\r
+    //\r
+    if (newLayer->fType==KMCLayer::kTypeNA) printf("Attention: the layer %s has undefined type\n",name);\r
+    //\r
+    newLayer->fIsDead =  (newLayer->fPhiRes==RIDICULOUS && newLayer->fZRes==RIDICULOUS);\r
+    //\r
+    if (fLayers.GetEntries()==0) \r
+      fLayers.Add(newLayer);\r
+    else {\r
+      //\r
+      for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
+       KMCLayer *l = (KMCLayer*)fLayers.At(i);\r
+       if (radius<l->fR) { fLayers.AddBefore(l,newLayer); break; }\r
+         if (radius>l->fR && (i+1)==fLayers.GetEntries() ) fLayers.Add(newLayer); // even bigger then last one\r
+      }\r
+      //\r
+    }\r
+    //\r
+    ClassifyLayers();\r
+    //\r
+  } else {\r
+    printf("Layer with the name %s does already exist\n",name);\r
+  }\r
+}\r
+\r
+//____________________________________________________________\r
+void KMCDetector::ClassifyLayers()\r
+{\r
+  // assign active Id's, etc\r
+  fLastActiveLayer = -1;\r
+  fLastActiveITSLayer = -1;\r
+  fNActiveLayers = 0;\r
+  fNActiveITSLayers = 0;\r
+  //\r
+  int nl = GetNLayers();\r
+  for (int il=0;il<nl;il++) {\r
+    KMCLayer* lr = GetLayer(il);\r
+    lr->SetUniqueID(il);\r
+    if (!lr->IsDead()) {\r
+      fLastActiveLayer = il; \r
+      lr->fActiveID = fNActiveLayers++;\r
+      if (lr->IsITS()) {\r
+       fLastActiveITSLayer = il;\r
+       fNActiveITSLayers++;\r
+      }\r
+    }\r
+  }\r
+  //\r
+  KMCProbe::SetNITSLayers(fNActiveITSLayers);\r
+}\r
+\r
+//____________________________________________________________\r
+void KMCDetector::KillLayer(char *name) {\r
+  //\r
+  // Marks layer as dead. Contribution only by Material Budget\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot mark as dead\n",name);\r
+  else {\r
+     tmp->fPhiRes = 999999;\r
+     tmp->fZRes = 999999;\r
+     tmp->fIsDead = kTRUE;\r
+     ClassifyLayers();\r
+  }\r
+}\r
+\r
+void KMCDetector::SetRadius(char *name, Float_t radius) {\r
+  //\r
+  // Set layer radius [cm]\r
+  //\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  //\r
+  if (!tmp) {\r
+    printf("Layer %s not found - cannot set radius\n",name);\r
+  } else {\r
+      \r
+    Float_t tmpRadL  = tmp->fx2X0;\r
+    Float_t tmpPhiRes = tmp->fPhiRes;\r
+    Float_t tmpZRes = tmp->fZRes;\r
+    Float_t tmpXRho = tmp->fXRho;\r
+    RemoveLayer(name); // so that the ordering is correct\r
+    AddLayer(name,radius,tmpRadL,tmpXRho,tmpPhiRes,tmpZRes);\r
+  }\r
+}\r
+\r
+Float_t KMCDetector::GetRadius(char *name) {\r
+  //\r
+  // Return layer radius [cm]\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot get radius\n",name);\r
+  else \r
+    return tmp->fR;\r
+\r
+  return 0;\r
+}\r
+\r
+void KMCDetector::SetRadiationLength(char *name, Float_t x2X0) {\r
+  //\r
+  // Set layer material [cm]\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot set layer material\n",name);\r
+  else {\r
+    tmp->fx2X0 = x2X0;\r
+  }\r
+}\r
+\r
+Float_t KMCDetector::GetRadiationLength(char *name) {\r
+  //\r
+  // Return layer radius [cm]\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot get layer material\n",name);\r
+  else \r
+    return tmp->fx2X0;\r
+    \r
+  return 0;\r
+  \r
+}\r
+\r
+void KMCDetector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {\r
+  //\r
+  // Set layer resolution in [cm]\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot set resolution\n",name);\r
+  else {\r
+    tmp->fPhiRes = phiRes;\r
+    tmp->fZRes = zRes;\r
+    tmp->fIsDead = (zRes==RIDICULOUS && phiRes==RIDICULOUS);\r
+    ClassifyLayers();\r
+  }\r
+}\r
+\r
+Float_t KMCDetector::GetResolution(char *name, Int_t axis) {\r
+  //\r
+  // Return layer resolution in [cm]\r
+  // axis = 0: resolution in rphi\r
+  // axis = 1: resolution in z\r
+  //\r
+\r
+  KMCLayer *tmp = GetLayer(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot get resolution\n",name);\r
+  else {\r
+    if (axis==0) return tmp->fPhiRes;\r
+    if (axis==1) return tmp->fZRes;\r
+    printf("error: axis must be either 0 or 1 (rphi or z axis)\n");\r
+  }\r
+  return 0;\r
+}\r
+\r
+void KMCDetector::SetLayerEfficiency(char *name, Float_t eff) {\r
+  //\r
+  // Set layer efficnecy (prop that his is missed within this layer) \r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot set layer efficiency\n",name);\r
+  else {\r
+    tmp->fEff = eff;\r
+  }\r
+}\r
+\r
+Float_t KMCDetector::GetLayerEfficiency(char *name) {\r
+  //\r
+  // Get layer efficnecy (prop that his is missed within this layer) \r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot get layer efficneicy\n",name);\r
+  else \r
+    return tmp->fEff;\r
+    \r
+  return 0;\r
+  \r
+}\r
+\r
+void KMCDetector::RemoveLayer(char *name) {\r
+  //\r
+  // Removes a layer from the list\r
+  //\r
+\r
+  KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
+  if (!tmp) \r
+    printf("Layer %s not found - cannot remove it\n",name);\r
+  else {\r
+    fLayers.Remove(tmp);\r
+    ClassifyLayers();\r
+  }\r
+}\r
+\r
+\r
+void KMCDetector::PrintLayout() {\r
+  //\r
+  // Prints the detector layout\r
+  //\r
+\r
+  printf("Detector %s: \"%s\"\n",GetName(),GetTitle());\r
+  \r
+  if (fLayers.GetEntries()>0) \r
+    printf("  Name \t\t r [cm] \t  X0 \t  phi & z res [um]\n");\r
+\r
+  KMCLayer *tmp = 0;\r
+  for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
+    tmp = (KMCLayer*)fLayers.At(i);\r
+  \r
+    // don't print all the tpc layers\r
+    TString name(tmp->GetName());\r
+    if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
+\r
+    printf("%d. %s \t %03.2f   \t%1.4f\t  ",i,\r
+          tmp->GetName(), tmp->fR, tmp->fx2X0);\r
+    if (tmp->fPhiRes==RIDICULOUS) \r
+      printf("  -  ");\r
+    else\r
+      printf("%3.0f   ",tmp->fPhiRes*10000);\r
+    if (tmp->fZRes==RIDICULOUS) \r
+      printf("  -\n");\r
+    else\r
+      printf("%3.0f\n",tmp->fZRes*10000);\r
+  }\r
+}\r
+\r
+void KMCDetector::PlotLayout(Int_t plotDead) {\r
+  //\r
+  // Plots the detector layout in Front view\r
+  //\r
+\r
+  Double_t x0=0, y0=0;\r
+\r
+  TGraphErrors *gr = new TGraphErrors();\r
+  gr->SetPoint(0,0,0);\r
+  KMCLayer *lastLayer = (KMCLayer*)fLayers.At(fLayers.GetEntries()-1);  Double_t maxRad = lastLayer->fR;\r
+  gr->SetPointError(0,maxRad,maxRad);\r
+  gr->Draw("APE");\r
+  \r
+\r
+  KMCLayer *tmp = 0;\r
+  for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {\r
+    tmp = (KMCLayer*)fLayers.At(i);\r
+  \r
+\r
+    Double_t txtpos = tmp->fR;\r
+    if ((tmp->IsDead())) txtpos*=-1; //\r
+    TText *txt = new TText(x0,txtpos,tmp->GetName());\r
+    txt->SetTextSizePixels(5); txt->SetTextAlign(21);\r
+    if (!tmp->IsDead() || plotDead) txt->Draw();\r
+\r
+    TEllipse *layEl = new TEllipse(x0,y0,tmp->fR);\r
+    //  layEl->SetFillColor(5);\r
+    layEl->SetFillStyle(5001);\r
+    layEl->SetLineStyle(tmp->IsDead()+1); // dashed if not active\r
+    layEl->SetLineColor(4);\r
+    TString name(tmp->GetName());\r
+    if (!tmp->IsDead()) layEl->SetLineWidth(2);\r
+    if (name.Contains("tpc") )  layEl->SetLineColor(29);\r
+\r
+    if (!tmp->IsDead() || plotDead) layEl->Draw();\r
+  \r
+  }\r
+\r
+}\r
+\r
+\r
+\r
+void KMCDetector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {\r
+  //\r
+  // Emulates the TPC\r
+  // \r
+  // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow \r
+\r
+\r
+  AddLayer((char*)"IFCtpc",   77.8,0.01367, 0); // Inner Field cage (RS: set correct xrho for eloss)\r
+  \r
+  // % Radiation Lengths ... Average per TPC row  (i.e. total/159 )\r
+  Float_t x2X0PerRow = 0.000036;\r
+  \r
+  Float_t tpcInnerRadialPitch  =    0.75 ;    // cm\r
+  Float_t tpcMiddleRadialPitch =    1.0  ;    // cm\r
+  Float_t tpcOuterRadialPitch  =    1.5  ;    // cm\r
+  //  Float_t tpcInnerPadWidth     =    0.4  ;    // cm\r
+  //  Float_t tpcMiddlePadWidth    =    0.6   ;   // cm\r
+  //  Float_t tpcOuterPadWidth     =    0.6   ;   // cm\r
+  Float_t innerRows            =   63 ;\r
+  Float_t middleRows           =   64  ;\r
+  Float_t outerRows            =   32  ;\r
+  Float_t tpcRows            =   (innerRows + middleRows + outerRows) ;\r
+  Float_t rowOneRadius         =   85.2  ;    // cm\r
+  Float_t row64Radius          =  135.1  ;    // cm\r
+  Float_t row128Radius         =  199.2  ;    // cm                       \r
\r
+  for ( Int_t k = 0 ; k < tpcRows ; k++ ) {\r
+    \r
+    Float_t rowRadius =0;\r
+    if (k<innerRows) \r
+      rowRadius =  rowOneRadius + k*tpcInnerRadialPitch ;\r
+    else if ( k>=innerRows && k<(innerRows+middleRows) )\r
+      rowRadius =  row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;\r
+    else if (k>=(innerRows+middleRows) && k<tpcRows )\r
+      rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;\r
+\r
+    if ( k%skip == 0 )\r
+      AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0, phiResMean,zResMean);    \r
+    else \r
+      AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0); // non "active" row\r
+    \r
+  \r
+  }\r
\r
+}\r
+\r
+void KMCDetector::RemoveTPC() {\r
+\r
+  // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)\r
+  KMCLayer *tmp = 0;\r
+  for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
+    tmp = (KMCLayer*)fLayers.At(i);  \r
+    TString name(tmp->GetName());\r
+    if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }\r
+  }\r
+  RemoveLayer((char*)"IFC");\r
+  \r
+}\r
+\r
+\r
+Double_t KMCDetector::ThetaMCS ( Double_t mass, Double_t x2X0, Double_t momentum ) const\r
+{\r
+  //\r
+  // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)\r
+  //\r
+\r
+  Double_t beta  =  momentum / TMath::Sqrt(momentum*momentum+mass*mass)  ;\r
+  Double_t theta =  0.0 ;    // Momentum and mass in GeV\r
+  // if ( RadLength > 0 ) theta  =  0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );\r
+  if ( x2X0 > 0 ) theta  =  0.0136 * TMath::Sqrt(x2X0) / ( beta * momentum ) * (1+0.038*TMath::Log(x2X0)) ;\r
+  return (theta) ;\r
+}\r
+\r
+\r
+Double_t KMCDetector::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
+{\r
+  // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm \r
+  // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm\r
+  // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius\r
+  Double_t sx, sy, goodHit ;\r
+  sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;\r
+  sy = 2 * TMath::Pi() *  searchRadiusZ    * searchRadiusZ    * HitDensity(radius) ;\r
+  goodHit =  TMath::Sqrt(1./((1+sx)*(1+sy)))  ;\r
+  return ( goodHit ) ;\r
+}\r
+\r
+\r
+Double_t KMCDetector::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
+{\r
+  // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm\r
+  // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function\r
+  Double_t sx, goodHit ;\r
+  sx = 2 * TMath::Pi() *  searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;\r
+  goodHit =  1./(1+sx) ;\r
+  return ( goodHit ) ;  \r
+}\r
+\r
+Double_t KMCDetector::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
+{\r
+  // Based on work by Ruben Shahoyen \r
+  // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function\r
+  // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account \r
+  // Following is correct for 2 DOF\r
+\r
+  Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
+  Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; \r
+  Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));\r
+  return ( goodHit ) ;  \r
+}\r
+\r
+Double_t KMCDetector::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
+{\r
+  // Based on work by Ruben Shahoyen \r
+  // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)\r
+\r
+  Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
+  Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; \r
+  Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));\r
+  return ( nullHit ) ;  \r
+}\r
+\r
+Double_t KMCDetector::HitDensity ( Double_t radius ) \r
+{\r
+  // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.\r
+  // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results\r
+  // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.\r
+  // Note that this function assumes we are working in CM and CM**2 [not meters].\r
+  // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.\r
+\r
+  //  Double_t MaxRadiusSlowDet = 0.1; //?   // Maximum radius for slow detectors.  Fast detectors \r
+                                        // and only fast detectors reside outside this radius.\r
+  Double_t arealDensity = 0 ;\r
+  if (radius<0.01) return 0;\r
+\r
+  if ( radius >= fMaxRadiusSlowDet ) \r
+    {\r
+      arealDensity  = OneEventHitDensity(fdNdEtaCent,radius)  ; // Fast detectors see central collision density (only)\r
+      arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius)  ;  // Increase density due to background \r
+    }\r
+\r
+  if (radius < fMaxRadiusSlowDet )\r
+    { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.\r
+      arealDensity  = OneEventHitDensity(fdNdEtaCent,radius) \r
+                   + IntegratedHitDensity(dNdEtaMinB,radius) \r
+                   + UpcHitDensity(radius) ;\r
+      arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;  \r
+      // Increase density due to background \r
+    } \r
+\r
+  return ( arealDensity ) ;  \r
+}\r
+\r
+\r
+double KMCDetector::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const\r
+{\r
+  // This is for one event at the vertex.  No smearing.\r
+  double den   = multiplicity / (2.*TMath::Pi()*radius*radius) * fDensFactorEta ; // 2 eta ?\r
+  // note: surface of sphere is  '4*pi*r^2'\r
+  //       surface of cylinder is '2*pi*r* h' \r
+  return den ;\r
+} \r
+\r
+\r
+double KMCDetector::IntegratedHitDensity(Double_t multiplicity, Double_t radius)\r
+{ \r
+  // The integral of minBias events smeared over a gaussian vertex distribution.\r
+  // Based on work by Yan Lu 12/20/2006, all radii in centimeters.\r
+\r
+  Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;\r
+  Double_t den   = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;\r
+\r
+  // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.\r
+  double dens1 = OneEventHitDensity(multiplicity,radius);\r
+  if ( den < dens1 )  den = dens1;\r
+\r
+  return den ;\r
+} \r
+\r
+\r
+double KMCDetector::UpcHitDensity(Double_t radius)\r
+{ \r
+  // QED electrons ...\r
+\r
+  Double_t mUPCelectrons ;                                 ;  \r
+  //  mUPCelectrons =  fLhcUPCscale * (1.23 - radius/6.5)      ;  // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC\r
+  mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB;      // Fit to 'Rossegger,Sadovsky'-Alice simulation\r
+  if ( mUPCelectrons < 0 ) mUPCelectrons =  0.0             ;  // UPC electrons fall off quickly and don't go to large R\r
+  mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ;  // UPCs increase Mulitiplicty ~ proportional to MinBias rate\r
+  mUPCelectrons *= UPCBackgroundMultiplier                 ;  // Allow for an external multiplier (eg 0-1) to turn off UPC\r
+\r
+  return mUPCelectrons ;\r
+} \r
+\r
+\r
+double KMCDetector::Dist(double z, double r)\r
+{\r
+  // Convolute dEta/dZ  distribution with assumed Gaussian of vertex z distribution\r
+  // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm\r
+  // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.\r
+  Int_t    index  =  1     ;     // Start weight at 1 for Simpsons rule integration\r
+  Int_t    nsteps =  301   ;     // NSteps must be odd for Simpson's rule to work\r
+  double   dist   =  0.0   ;\r
+  double   dz0    =  ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1)  ;  //cm\r
+  double    z0    =  0.0   ;     //cm\r
+  for(int i=0; i<nsteps; i++){\r
+    if ( i == nsteps-1 ) index = 1 ;\r
+    z0 = -4*SigmaD + i*dz0 ;\r
+    dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) * \r
+      (1/sqrt((z-z0)*(z-z0) + r*r)) ;\r
+    if ( index != 4 ) index = 4; else index = 2 ;\r
+  }\r
+  return dist; \r
+}\r
+\r
+#define  PZero   0.861  // Momentum of back to back decay particles in the CM frame\r
+#define  EPiZero 0.872  // Energy of the pion from a D0 decay at rest\r
+#define  EKZero  0.993  // Energy of the Kaon from a D0 decay at rest\r
+\r
+Double_t KMCDetector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) const {\r
+  // Math from Ron Longacre.  Note hardwired energy to bin conversion for PtK and PtPi.\r
+\r
+  Double_t const1  =  pt / D0Mass ;\r
+  Double_t const2  =  TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;\r
+  Double_t sum, ptPi, ptK ;\r
+  Double_t effp, effk ;\r
+\r
+  sum = 0.0 ;\r
+  for ( Int_t k = 0 ; k < 360 ; k++ )   {\r
+    \r
+    Double_t theta = k * TMath::Pi() / 180. ;\r
+    \r
+    ptPi = TMath::Sqrt( \r
+                      PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +\r
+                      const1*const1*EPiZero*EPiZero -\r
+                      2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +\r
+                      PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)\r
+                      ) ;\r
+    \r
+    ptK = TMath::Sqrt( \r
+                     PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +\r
+                     const1*const1*EKZero*EKZero +\r
+                     2*PZero*TMath::Cos(theta)*const2*const1*EKZero +\r
+                     PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)\r
+                     ) ;\r
+\r
+    // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays\r
+    Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ; \r
+    Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ; \r
+      \r
+    if ( pionindex >= 20 ) pionindex = 399 ;\r
+    if ( pionindex >= 0 )   effp = corrEfficiency[0][pionindex] ;\r
+    if ( pionindex <  0 )   effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd\r
+    if ( effp < 0 )         effp = 0 ;\r
+\r
+    if ( kaonindex >= 20 ) kaonindex = 399 ;\r
+    if ( kaonindex >= 0 )   effk = corrEfficiency[1][kaonindex] ;\r
+    if ( kaonindex <  0 )   effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd\r
+    if ( effk < 0 )         effk = 0 ;\r
+\r
+    // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.\r
+      \r
+    sum += effp * effk ;\r
\r
+  }    \r
+  \r
+  Double_t mean =sum/360; \r
+  return mean ;\r
+  \r
+}\r
+\r
+KMCProbe* KMCDetector::PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi, double x,double y, double z)\r
+{\r
+  // Prepare trackable Kalman track at the farthest position\r
+  //\r
+  // Set track parameters\r
+  // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis\r
+  fProbe.Reset();\r
+  fProbe.SetMass(mass);\r
+  KMCProbe* probe = new KMCProbe(fProbe);\r
+  double *trPars = (double*)probe->GetParameter();\r
+  double *trCov  = (double*)probe->GetCovariance();\r
+  double xyz[3] = {x,y,z};\r
+  probe->Global2LocalPosition(xyz,phi);\r
+  probe->Set(xyz[0],phi,trPars,trCov);\r
+  trPars[KMCProbe::kY] = xyz[1];\r
+  trPars[KMCProbe::kZ] = xyz[2];\r
+  trPars[KMCProbe::kSnp] = 0;                       //            track along X axis at the vertex\r
+  trPars[KMCProbe::kTgl] = TMath::Tan(lambda);                // dip\r
+  trPars[KMCProbe::kPtI] = charge/pt;               //            q/pt      \r
+  //\r
+  // put tiny errors to propagate to the outer-most radius\r
+  trCov[KMCProbe::kY2] = trCov[KMCProbe::kZ2] = trCov[KMCProbe::kSnp2] = trCov[KMCProbe::kTgl2] = trCov[KMCProbe::kPtI2] = 1e-20;\r
+  fProbe = *probe;  // store original track\r
+  //\r
+  // propagate to last layer\r
+  fLastActiveLayerTracked = 0;\r
+  for (Int_t j=0; j<=fLastActiveLayer; j++) {\r
+    KMCLayer* lr = GetLayer(j);\r
+    lr->Reset();\r
+    //\r
+    if (!PropagateToLayer(probe,lr,1)) break;\r
+    if (!probe->CorrectForMeanMaterial(lr, kFALSE)) break;\r
+    //\r
+    lr->fClCorr.Set(probe->GetY(),probe->GetZ(), probe->GetX(), probe->GetAlpha());\r
+    if (!lr->IsDead()) fLastActiveLayerTracked = j;\r
+  }\r
+  probe->ResetCovMat();// reset cov.matrix\r
+  printf("Last active layer trracked: %d (out of %d)\n",fLastActiveLayerTracked,fLastActiveLayer);\r
+  //\r
+  return probe;\r
+}\r
+\r
+\r
+\r
+TGraph * KMCDetector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {\r
+  //\r
+  // returns the momentum resolution \r
+  //\r
+  \r
+  TGraph *graph = new TGraph(20, fTransMomenta, fMomentumRes);\r
+  graph->SetTitle("Momentum Resolution .vs. Pt" ) ;\r
+  //  graph->GetXaxis()->SetRangeUser(0.,5.0) ;\r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;\r
+  graph->GetYaxis()->CenterTitle();\r
+\r
+  graph->SetMaximum(20) ;\r
+  graph->SetMinimum(0.1) ;\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+\r
+}\r
+\r
+TGraph * KMCDetector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {\r
\r
+  // Returns the pointing resolution\r
+  // axis = 0 ... rphi pointing resolution\r
+  // axis = 1 ... z pointing resolution\r
+  //\r
+\r
+  TGraph * graph =  0;\r
+\r
+  if (axis==0) {\r
+    graph = new TGraph ( 20, fTransMomenta, fResolutionRPhi ) ;\r
+    graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;\r
+    graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;\r
+  } else {\r
+    graph =  new TGraph ( 20, fTransMomenta, fResolutionZ ) ;\r
+    graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;\r
+    graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;\r
+  }\r
+  \r
+  graph->SetMinimum(1) ;\r
+  graph->SetMaximum(300.1) ;\r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->CenterTitle();\r
+  \r
+  graph->SetLineWidth(linewidth);\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  \r
+  return graph;\r
+\r
+}\r
+\r
+\r
+TGraph * KMCDetector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {\r
+  //\r
+  // returns the Pointing resolution (accoring to Telescope equation)\r
+  // axis =0 ... in rphi\r
+  // axis =1 ... in z\r
+  //\r
+  \r
+  Double_t resolution[20];\r
+\r
+  Double_t layerResolution[2];\r
+  Double_t layerRadius[2];\r
+  Double_t layerThickness[2];\r
+\r
+  Int_t count =0; // search two first active layers\r
+  printf("Telescope equation for layers:  ");\r
+  for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
+    KMCLayer *l = (KMCLayer*)fLayers.At(i);\r
+    if (!l->IsDead() && l->fR>0) {\r
+      layerRadius[count]     = l->fR;\r
+      layerThickness[count]  = l->fx2X0;\r
+      if (axis==0) {\r
+       layerResolution[count] = l->fPhiRes;\r
+      } else {\r
+       layerResolution[count] = l->fZRes;\r
+      }\r
+      printf("%s, ",l->GetName());\r
+      count++;\r
+    }\r
+    if (count>=2) break;       \r
+  }\r
+  printf("\n");\r
+\r
+  Double_t pt, momentum, thickness,aMCS ;\r
+  Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
+\r
+  for ( Int_t i = 0 ; i < 20 ; i++ ) { \r
+    // Reference data as if first two layers were acting all alone \r
+    pt  =  fTransMomenta[i]  ;\r
+    momentum = pt / TMath::Cos(lambda)   ;  // Total momentum\r
+    resolution[i] =  layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1] \r
+      +  layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;\r
+    resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;\r
+    thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;\r
+    aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;\r
+    resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;\r
+    resolution[i] =  TMath::Sqrt(resolution[i]) * 10000.0 ;  // result in microns\r
+  }\r
+\r
+\r
+\r
+  TGraph* graph = new TGraph ( 20, fTransMomenta, resolution ) ;\r
+   \r
+  if (axis==0) {\r
+    graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;\r
+    graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;\r
+  } else {\r
+    graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;\r
+    graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;\r
+  }\r
+  graph->SetMinimum(1) ;\r
+  graph->SetMaximum(300.1) ;\r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->CenterTitle();\r
+  \r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineStyle(kDashed);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+\r
+}\r
+\r
+TGraph * KMCDetector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {\r
+  //\r
+  // particle = 0 ... choosen particle (setted particleMass)\r
+  // particle = 1 ... Pion\r
+  // particle = 2 ... Kaon\r
+  // particle = 3 ... D0\r
+  //\r
+  Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
+  \r
+  Double_t particleEfficiency[20]; // with chosen particle mass\r
+  Double_t kaonEfficiency[20], pionEfficiency[20], d0efficiency[20]; \r
+  Double_t partEfficiency[2][20];\r
+  \r
+  if (particle != 0) {\r
+    // resulting Pion and Kaon efficiency scaled with overall efficiency\r
+    Double_t doNotDecayFactor;\r
+    for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon\r
+      \r
+      for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+       // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.\r
+       Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity\r
+       if ( massloop == 1 ) { // KAON\r
+         doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.\r
+         kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;\r
+       } else { // PION\r
+         doNotDecayFactor = 1.0 ;\r
+         pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;       \r
+       }\r
+       partEfficiency[0][j] = pionEfficiency[j];\r
+       partEfficiency[1][j] = kaonEfficiency[j];\r
+      }      \r
+    }\r
+    \r
+    // resulting estimate of the D0 efficiency\r
+    for ( Int_t j = 0 ; j < 20 ; j++ ) {\r
+      d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);\r
+    }\r
+  } else { \r
+    for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+      particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;\r
+      // NOTE: Decay factor (see kaon) should be included to be realiable\r
+    }\r
+  }\r
+\r
+  for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+    pionEfficiency[j]     *= 100;\r
+    kaonEfficiency[j]     *= 100;\r
+    d0efficiency[j]       *= 100;\r
+    particleEfficiency[j] *= 100;\r
+  }\r
\r
+  TGraph * graph =  0;\r
+  if (particle==0) {\r
+    graph = new TGraph ( 20, fTransMomenta, particleEfficiency ) ; // choosen mass\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle==1) {\r
+    graph = new TGraph ( 20, fTransMomenta, pionEfficiency ) ;\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle ==2) {\r
+    graph = new TGraph ( 20, fTransMomenta, kaonEfficiency ) ;\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle ==3) {\r
+    graph = new TGraph ( 20, fTransMomenta, d0efficiency ) ;\r
+    graph->SetLineStyle(kDashed);\r
+  } else \r
+    return 0;\r
+\r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->SetTitle("Efficiency (%)") ;\r
+  graph->GetYaxis()->CenterTitle();\r
+         \r
+  graph->SetMinimum(0.01) ; \r
+  graph->SetMaximum(100)  ; \r
+\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+}\r
+\r
+TGraph * KMCDetector::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {\r
+  //\r
+  // particle = 0 ... choosen particle (setted particleMass)\r
+  // particle = 1 ... Pion\r
+  // particle = 2 ... Kaon\r
+  //\r
+\r
+  Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
+  \r
+  Double_t particleFake[20]; // with chosen particle mass\r
+  Double_t kaonFake[20], pionFake[20];\r
+  Double_t partFake[2][20];\r
+  \r
+  if (particle != 0) {\r
+    // resulting Pion and Kaon efficiency scaled with overall efficiency\r
+    Double_t doNotDecayFactor;\r
+    for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon\r
+      \r
+      for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+       // JT Test Let the kaon decay.  If it decays inside the TPC ... then it is gone; for all decays < 130 cm.\r
+       Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda)           ;  // Total momentum at average rapidity\r
+       if ( massloop == 1 ) { // KAON\r
+         doNotDecayFactor  = TMath::Exp(-130/(371*momentum/KaonMass)) ;  // Decay length for kaon is 371 cm.\r
+         kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;\r
+       } else { // PION\r
+         pionFake[j] = fFake[0][j] ;   \r
+       }\r
+       partFake[0][j] = pionFake[j];\r
+       partFake[1][j] = kaonFake[j];\r
+      }      \r
+    }\r
+    \r
+  } else { \r
+    for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+      particleFake[j] = fFake[2][j];\r
+      // NOTE: Decay factor (see kaon) should be included to be realiable\r
+    }\r
+  }\r
+\r
+  for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
+    pionFake[j]     *= 100;\r
+    kaonFake[j]     *= 100;\r
+    particleFake[j] *= 100;\r
+  }\r
\r
+  TGraph * graph =  0;\r
+  if (particle==0) {\r
+    graph = new TGraph ( 20, fTransMomenta, particleFake ) ; // choosen mass\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle==1) {\r
+    graph = new TGraph ( 20, fTransMomenta, pionFake ) ;\r
+    graph->SetLineWidth(1);\r
+  }  else if (particle ==2) {\r
+    graph = new TGraph ( 20, fTransMomenta, kaonFake ) ;\r
+    graph->SetLineWidth(1);\r
+  } \r
+  \r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->SetTitle("Fake (%)") ;\r
+  graph->GetYaxis()->CenterTitle();\r
+         \r
+  graph->SetMinimum(0.01) ; \r
+  graph->SetMaximum(100)  ; \r
+\r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+}\r
+\r
+\r
+TGraph* KMCDetector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {\r
+  //\r
+  // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)\r
+  // mode 0: impact parameter (convolution of pointing and vertex resolution)\r
+  // mode 1: pointing resolution\r
+  // mode 2: vtx resolution \r
+  \r
+  \r
+  TGraph *graph = new TGraph();\r
+\r
+  //  TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?\r
+  TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); // \r
+  TFormula vtxResZ("vtxResZ","600/(x+6)+10"); // \r
+    \r
+  TGraph *trackRes = GetGraphPointingResolution(axis,1);\r
+  Double_t *pt = trackRes->GetX();\r
+  Double_t *trRes = trackRes->GetY();\r
+  for (Int_t ip =0; ip<trackRes->GetN(); ip++) {\r
+    Double_t vtxRes = 0;\r
+    if (axis==0) \r
+      vtxRes = vtxResRPhi.Eval(pt[ip]);\r
+    else \r
+      vtxRes = vtxResZ.Eval(pt[ip]);\r
+    \r
+    if (mode==0)\r
+      graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));\r
+    else if (mode ==1)\r
+      graph->SetPoint(ip,pt[ip],trRes[ip]);\r
+    else\r
+      graph->SetPoint(ip,pt[ip],vtxRes);\r
+  }\r
+  \r
+  graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;\r
+  graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;\r
+  \r
+  graph->SetMinimum(1) ;\r
+  graph->SetMaximum(300.1) ;\r
+  graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
+  graph->GetXaxis()->CenterTitle();\r
+  graph->GetXaxis()->SetNoExponent(1) ;\r
+  graph->GetXaxis()->SetMoreLogLabels(1) ;\r
+  graph->GetYaxis()->CenterTitle();\r
+  \r
+  graph->SetLineColor(color);\r
+  graph->SetMarkerColor(color);\r
+  graph->SetLineWidth(linewidth);\r
+\r
+  return graph;\r
+\r
+}\r
+\r
+TGraph* KMCDetector::GetGraph(Int_t number, Int_t color, Int_t linewidth) {\r
+  // \r
+  // returns graph according to the number\r
+  //\r
+  switch(number) {\r
+  case 1:\r
+    return GetGraphPointingResolution(0,color, linewidth); // dr\r
+  case 2:\r
+    return GetGraphPointingResolution(1,color, linewidth); // dz\r
+  case 3:\r
+    return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele\r
+  case 4:\r
+    return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele\r
+  case 5:\r
+    return GetGraphMomentumResolution(color, linewidth); // pt resolution\r
+  case 10:\r
+    return GetGraphRecoEfficiency(0, color, linewidth);  // tracked particle\r
+  case 11:\r
+    return GetGraphRecoEfficiency(1, color, linewidth);  // eff. pion\r
+  case 12:\r
+    return GetGraphRecoEfficiency(2, color, linewidth);  // eff. kaon\r
+  case 13: \r
+    return GetGraphRecoEfficiency(3, color, linewidth);  // eff. D0\r
+  case 15:\r
+    return GetGraphRecoFakes(0, color, linewidth);  // Fake tracked particle\r
+  case 16:\r
+    return GetGraphRecoFakes(1, color, linewidth);  // Fake pion\r
+  case 17:\r
+    return GetGraphRecoFakes(2, color, linewidth);  // Fake kaon\r
+  default:\r
+    printf(" Error: chosen graph number not valid\n");\r
+  }\r
+  return 0;\r
+\r
+}\r
+\r
+void KMCDetector::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon, int setVer) {\r
+  \r
+  // All New configuration with X0 = 0.3 and resolution = 4 microns\r
+  \r
+  AddLayer((char*)"bpipe_its",2.0,0.0022, 0.092); // beam pipe, 0.5 mm Be\r
+  AddLayer((char*)"vertex_its",     0,     0); // dummy vertex for matrix calculation\r
+  if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {\r
+    printf("vertex will work as constraint: %.4f %.4f\n",fgVtxConstraint[0],fgVtxConstraint[1]);\r
+    SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);\r
+  }\r
+  //\r
+  // new ideal Pixel properties?\r
+  Double_t x0     = 0.0050;\r
+  Double_t resRPhi = 0.0006;\r
+  Double_t resZ   = 0.0006;\r
+  Double_t xrho = 0.0116;  // assume 0.5mm of silicon for eloss\r
+  //\r
+  if (flagMon) {\r
+    x0 *= 3./5.;\r
+    xrho *=3./5.;\r
+    resRPhi = 0.0004;\r
+    resZ   = 0.0004;\r
+  }\r
+\r
+  double sclD = 1;\r
+  double sclZ = 1;\r
+\r
+  // default all new  \r
+  if (setVer<=0) {\r
+    AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ); \r
+  }\r
+  else if (setVer==1) {\r
+    AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its2",  2.8 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its3",  3.6 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its4", 20.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its5", 22.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its6", 43.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its7", 43.6 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    //\r
+    /*\r
+    UInt_t patt[] = {\r
+      KMCTrackSummary::Bits(1,1,1),\r
+      KMCTrackSummary::Bits(0,0,0,1,1),\r
+      KMCTrackSummary::Bits(0,0,0,0,0,1,1)\r
+    };\r
+    RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\r
+    */\r
+  }\r
+  else if (setVer==2) {\r
+    AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its1a", 2.8 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its2",  3.6 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its2a", 4.2 ,  x0, xrho, resRPhi, resZ); \r
+    AddLayer((char*)"its3", 20.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its4", 22.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its5", 33.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its6", 43.0 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    AddLayer((char*)"its7", 43.6 ,  x0, xrho, resRPhi*sclD, resZ*sclZ); \r
+    //\r
+    /*\r
+    UInt_t patt[] = {\r
+      KMCTrackSummary::Bits(1,1,1,1),\r
+      KMCTrackSummary::Bits(0,0,0,0,1,1),\r
+      KMCTrackSummary::Bits(0,0,0,0,0,0,1,1,1)\r
+    };\r
+    RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\r
+    */\r
+  }\r
+   /*\r
+  // last 2 layers strips\r
+  double resRPStr=0.0020, resZStr = 0.0830;\r
+  AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPStr, resZStr); \r
+  AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPStr, resZStr); \r
+   */\r
+  //*/\r
+  /*\r
+  // resolution scaled with R as res(2.2) * R/2.2\r
+  AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ); \r
+  KMCLayer* lr0 = 0;\r
+  for (int i=0;i<GetNActiveITSLayers();i++) {\r
+    KMCLayer *lr = GetActiveLayer(i);\r
+    if (lr->IsVertex()) continue;\r
+    if (lr0==0) {\r
+      lr0=lr; \r
+      //  continue;\r
+    }\r
+    double scl = 5*lr->GetRadius()/lr0->GetRadius();\r
+    SetResolution((char*)lr->GetName(), resRPhi*scl, resZ*scl*4);\r
+  }\r
+  */\r
+  /*\r
+  // 1st 2 layers double\r
+  AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its1a",  2.8 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its2a",  4.2 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its3a",  6.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ); \r
+  */\r
+  /*\r
+  // last 4 layers doubled\r
+  AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its4a", 12.8 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its5", 23.5 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its5a", 23.9 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its6", 39.6 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its6a", 40.0 ,  x0, xrho, resRPhi, resZ); \r
+\r
+  AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its7a", 43.4 ,  x0, xrho, resRPhi, resZ); \r
+  */\r
+\r
+  /* //last 3 lr together \r
+  AddLayer((char*)"its1",  2.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its2",  3.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its3",  6.8 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its4", 12.4 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its5", 42.2 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its6", 42.6 ,  x0, xrho, resRPhi, resZ); \r
+  AddLayer((char*)"its7", 43.0 ,  x0, xrho, resRPhi, resZ); \r
+  */\r
+  if (flagTPC) {\r
+    AddTPC(0.1,0.1);                        // TPC\r
+  }\r
+  PrintLayout();\r
+}\r
+\r
+void KMCDetector::MakeAliceCurrent(Bool_t flagTPC, Int_t AlignResiduals) {\r
+\r
+  // Numbers taken from \r
+  // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks\r
+  // number for misalingment: private communication with Andrea Dainese\r
+\r
+  //  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};\r
+\r
+  AddLayer((char*)"vertex_its",     0,     0); // dummy vertex for matrix calculation\r
+  if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {\r
+    printf("vertex will work as constraint: %.4f %.4f",fgVtxConstraint[0],fgVtxConstraint[1]);\r
+    SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);\r
+  }\r
+  //\r
+  AddLayer((char*)"bpipe_its",2.94,0.0022, 1.48e-01); // beam pipe\r
+  AddLayer((char*)"tshld1_its",11.5,0.0065, 1.34e-01); // Thermal shield  // 1.3% /2\r
+  AddLayer((char*)"tshld2_its",31.0,0.0108, 2.22e-01); // Thermal shield  // 1.3% /2\r
+  //\r
+  if (flagTPC) {\r
+    AddTPC(0.1,0.1);                        // TPC\r
+  }\r
+  // Adding the ITS - current configuration\r
+  \r
+  if (AlignResiduals==0) {\r
+    AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, 0.0012, 0.0130);\r
+    AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, 0.0012, 0.0130);\r
+    AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, 0.0035, 0.0025);\r
+    AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, 0.0035, 0.0025);\r
+    AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 0.0020, 0.0830);\r
+    AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 0.0020, 0.0830);\r
+    /*\r
+    UInt_t patt[] = {\r
+      KMCTrackSummary::Bits(1,1),\r
+      KMCTrackSummary::Bits(0,0,1,1),\r
+      KMCTrackSummary::Bits(0,0,0,0,1,1)\r
+    };\r
+    RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\r
+    */\r
+  } else if (AlignResiduals==1) {\r
+\r
+    // tracking errors ...\r
+    // (Additional systematic errors due to misalignments) ... \r
+    // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]\r
+    // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);\r
+\r
+    AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010), \r
+            TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
+    AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),\r
+            TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
+    AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),\r
+            TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
+    AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),\r
+            TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
+    AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), \r
+            TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));\r
+    AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),\r
+            TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));   \r
+    \r
+  } else if (AlignResiduals==2) {\r
+\r
+    \r
+    // tracking errors ... PLUS ... module misalignment\r
+    \r
+    // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020);  // [cm]\r
+    // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);\r
+    \r
+    //  the ITS modules are misalignment with small gaussian smearings with\r
+    //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD\r
+    \r
+    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), \r
+            TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
+    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),\r
+            TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
+    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),\r
+            TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
+    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),\r
+            TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
+    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), \r
+            TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));\r
+    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),\r
+            TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); \r
+\r
+  } else {\r
+      \r
+      //  the ITS modules are misalignment with small gaussian smearings with\r
+      //  sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD\r
+      //  unknown in Z ????\r
+\r
+    AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), \r
+            TMath::Sqrt(0.0130*0.0130+0.000*0.000));\r
+    AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),\r
+            TMath::Sqrt(0.0130*0.0130+0.000*0.000));\r
+    AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),\r
+            TMath::Sqrt(0.0025*0.0025+0.000*0.000));\r
+    AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),\r
+            TMath::Sqrt(0.0025*0.0025+0.000*0.000));\r
+    AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), \r
+            TMath::Sqrt(0.0830*0.0830+0.000*0.000));\r
+    AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),\r
+            TMath::Sqrt(0.0830*0.0830+0.000*0.000));       \r
+  }\r
+  \r
+}\r
+\r
+\r
+void KMCDetector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {\r
+  //\r
+  // Produces the standard performace plots\r
+  //\r
\r
+  if (!add) {\r
+\r
+    TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);  \r
+    c1->Divide(2,2);\r
+    \r
+    c1->cd(1);  gPad->SetGridx();   gPad->SetGridy(); \r
+    gPad->SetLogx(); \r
+    TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);\r
+    eff->SetTitle("Efficiencies");\r
+    eff->Draw("AL");\r
+    if (!onlyPionEff) {\r
+      GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");\r
+      GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");\r
+    }\r
+    c1->cd(2); gPad->SetGridx();   gPad->SetGridy(); \r
+    gPad->SetLogy();  gPad->SetLogx(); \r
+    GetGraphMomentumResolution(color,linewidth)->Draw("AL");\r
+    \r
+    c1->cd(3); gPad->SetGridx();   gPad->SetGridy(); \r
+    gPad->SetLogx(); \r
+    GetGraphPointingResolution(0,color,linewidth)->Draw("AL");\r
+    \r
+    c1->cd(4); gPad->SetGridx();   gPad->SetGridy(); \r
+    gPad->SetLogx(); \r
+    GetGraphPointingResolution(1,color,linewidth)->Draw("AL");\r
+\r
+  } else {\r
+\r
+    TVirtualPad *c1 = gPad->GetMother();\r
+\r
+    c1->cd(1);\r
+    GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");\r
+    if (!onlyPionEff) {\r
+      GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");\r
+      GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");\r
+    }\r
+    c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");\r
+    \r
+    c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");\r
+    \r
+    c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");\r
+    \r
+  }\r
+\r
+}\r
+\r
+void KMCDetector::ApplyMS(KMCProbe* trc, double x2X0) const\r
+{\r
+  // simulate random modification of track params due to the MS\r
+  if (x2X0<=0) return;\r
+  double alpha = trc->GetAlpha(); // store original alpha\r
+  double mass = trc->GetMass();\r
+  //\r
+  double snp = trc->GetSnp();\r
+  double dip = trc->GetTgl();\r
+  Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));\r
+  x2X0 *= angle;\r
+  //\r
+  static double covCorr[15],covDum[21]={0};\r
+  static double mom[3],pos[3];\r
+  double *cov = (double*) trc->GetCovariance();\r
+  memcpy(covCorr,cov,15*sizeof(double));\r
+  trc->GetXYZ(pos);\r
+  trc->GetPxPyPz(mom);\r
+  double pt2 = mom[0]*mom[0]+mom[1]*mom[1];\r
+  double pt = TMath::Sqrt(pt2);\r
+  double ptot2 = pt2 + mom[2]*mom[2];\r
+  double ptot  = TMath::Sqrt(ptot2);\r
+  double beta = ptot/TMath::Sqrt(ptot2 + mass*mass);\r
+  double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);\r
+  //\r
+  // a la geant\r
+  double phiSC = gRandom->Rndm()*TMath::Pi();\r
+  double thtSC = gRandom->Gaus(0,1.4142*sigth);\r
+  //  printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);\r
+  double sn = TMath::Sin(thtSC);\r
+  double dx = sn*TMath::Sin(phiSC);\r
+  double dy = sn*TMath::Cos(phiSC);  \r
+  double dz = TMath::Cos(thtSC);\r
+  double v[3];\r
+  //  printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);\r
+  for (int i=3;i--;) mom[i] /= ptot;\r
+  double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);\r
+  if (!IsZero(pt)) {\r
+    double pd1 = mom[0]/vmm;\r
+    double pd2 = mom[1]/vmm;\r
+    v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;\r
+    v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;\r
+    v[2] = -vmm*dx                + mom[2]*dz;\r
+  }\r
+  else {\r
+    v[0] = dx;\r
+    v[1] = dy;\r
+    v[2] = dz*TMath::Sign(1.,mom[2]);\r
+  }\r
+  double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);\r
+  //  printf("before :%+e %+e %+e  || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2],  sigth, x2X0, pt, beta);\r
+  //  trc->Print();\r
+  // direction cosines -> p\r
+  for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;\r
+  //  printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);\r
+  trc->Set(pos,mom,covDum,trc->Charge());\r
+  //\r
+  trc->Rotate(alpha);\r
+  memcpy(cov,covCorr,15*sizeof(double));\r
+  //\r
+}\r
+\r
+//________________________________________________________________________________\r
+Bool_t KMCDetector::TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const\r
+{\r
+  // Transport track till layer maxLr, applying random MS\r
+  //\r
+  for (Int_t j=0; j<maxLr; j++) {\r
+    KMCLayer* lr0 = (KMCLayer*)fLayers.At(j);\r
+    KMCLayer* lr  = (KMCLayer*)fLayers.At(j+1);\r
+    if (!lr0->IsVertex()) {\r
+      ApplyMS(probTr,lr0->fx2X0); // apply MS\r
+      if (!probTr->CorrectForMeanMaterial(lr0,kFALSE)) return kFALSE; \r
+    }\r
+    //\r
+    if (!PropagateToLayer(probTr,lr,1)) return kFALSE;\r
+    // store randomized cluster local coordinates and phi\r
+    lr->ResetBgClusters();\r
+    double rz,ry;\r
+    gRandom->Rannor(rz,ry);\r
+    lr->GetMCCluster()->Set(probTr->GetY()+ry*lr->GetPhiRes(),probTr->GetZ()+rz*lr->GetZRes(), \r
+                              probTr->GetX(), probTr->GetAlpha() );    \r
+    //\r
+  }\r
+  //\r
+  return kTRUE;\r
+}\r
+\r
+//________________________________________________________________________________\r
+Bool_t KMCDetector::SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr,int nMC, int offset)\r
+{\r
+  // analityc and fullMC (nMC trials) evaluaion of tracks with given kinematics.\r
+  // the results are filled in KMCTrackSummary objects provided via summArr array\r
+  //\r
+  int progressP = 10;//0; // print progress in percents\r
+  //\r
+  progressP = int(nMC*0.01*progressP);\r
+\r
+  if (!SolveSingleTrackViaKalman(mass,pt,eta)) return kFALSE;\r
+  //\r
+  // Store non-updated track errors of inward propagated seed >>>>>>>>\r
+  int maxLr = fLastActiveITSLayer + offset;\r
+  if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;\r
+  KMCProbe probeTmp = fProbe; // original probe at vertex\r
+  KMCLayer* lr = 0;\r
+  for (Int_t j=1; j<=maxLr; j++) {\r
+    lr = GetLayer(j);\r
+    //    printf("Here0: %d\n",j);\r
+    if (!PropagateToLayer(&probeTmp,lr,1)) return 0;\r
+    if (j!=maxLr) if (!probeTmp.CorrectForMeanMaterial(lr, kFALSE)) return 0;\r
+    //    printf("Prelim. Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(probeTmp.GetSigmaY2()),TMath::Sqrt(probeTmp.GetSigmaZ2()));\r
+  }\r
+  for (Int_t j=maxLr; j>0; j--) {\r
+    lr = GetLayer(j);\r
+    //    printf("Here1: %d\n",j);\r
+    if (j!=maxLr) if (!PropagateToLayer(&probeTmp,lr,-1)) return 0;\r
+    lr->fSig2EstD = probeTmp.GetSigmaY2();\r
+    lr->fSig2EstZ = probeTmp.GetSigmaZ2();\r
+    //    probeTmp.Print("l");\r
+    printf("Natural Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(lr->fSig2EstD),TMath::Sqrt(lr->fSig2EstZ));\r
+    if (!probeTmp.CorrectForMeanMaterial(lr, kTRUE)) return 0;\r
+  }\r
+  // Store non-updated track errors of inward propagated seed <<<<<<<<\r
+  //\r
+  int nsm = sumArr ? sumArr->GetEntriesFast() : 0;\r
+  KMCLayer* vtx = GetLayer(0);\r
+  //\r
+  for (int i=0;i<nsm;i++) {\r
+    KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(i);\r
+    if (!tsm) continue;\r
+    tsm->SetRefProbe( GetProbeTrack() ); // attach reference track (generated)\r
+    tsm->SetAnProbe( vtx->GetAnProbe() ); // attach analitycal solution\r
+  }\r
+  //\r
+  TStopwatch sw;\r
+  sw.Start();\r
+  for (int it=0;it<nMC;it++) {\r
+    //    printf("ev: %d\n",it);\r
+    SolveSingleTrackViaKalmanMC(offset);\r
+    KMCProbe* trc = vtx->GetWinnerMCTrack();\r
+    if (progressP==1 || (progressP>0 &&  (it%progressP)==0)) {\r
+      printf("%d%% done |",it*100/nMC); \r
+      sw.Stop(); sw.Print(); sw.Start(kFALSE);\r
+    }\r
+    for (int ism=nsm;ism--;) { // account the track in each of summaries\r
+      KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(ism);\r
+      if (!tsm) continue;\r
+      tsm->AddUpdCalls(GetUpdCalls());\r
+      tsm->AddTrack(trc); \r
+    }\r
+  }\r
+  //\r
+  sw.Stop();\r
+  printf("Total time: "); sw.Print();\r
+  return kTRUE;\r
+}\r
+\r
+//________________________________________________________________________________\r
+Int_t KMCDetector::GetLayerID(int actID) const\r
+{\r
+  // find physical layer id from active id\r
+  if (actID<0 || actID>fNActiveLayers) return -1;\r
+  int start = actID<fNActiveITSLayers ? fLastActiveITSLayer : fLastActiveLayer;\r
+  for (int i=start+1;i--;) {\r
+    if (GetLayer(i)->GetActiveID()==actID) return i;   \r
+  }\r
+  return -1;\r
+}\r
+\r
+//________________________________________________________________________________\r
+KMCProbe* KMCDetector::KalmanSmooth(int actLr, int actMin,int actMax) const\r
+{\r
+  // estimate kalman smoothed track params at given active lr \r
+  // from fit at layers actMin:actMax (excluding actLr)\r
+  // SolveSingleTrackViaKalman must have been called before\r
+  //\r
+  if (actMin>actMax) swap(actMin,actMax);\r
+  if (actMax>=fNActiveLayers) actMax = fNActiveLayers-1;\r
+  int nlrfit = actMax-actMin;\r
+  if (actLr>=actMin && actLr<=actMax) nlrfit-=1;\r
+  if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}\r
+  static KMCProbe iwd,owd;\r
+  //\r
+  // find phisical layer id's\r
+  int pLr  = GetLayerID(actLr);\r
+  int pMin = GetLayerID(actMin);\r
+  int pMax = GetLayerID(actMax);\r
+  //\r
+  //  printf(">>> %d %d %d\n",pLr, pMin,pMax);\r
+  Bool_t useIwd=kFALSE, useOwd=kFALSE;\r
+  if (pLr<pMax) { // need inward piece\r
+    iwd = GetLayer(pMax)->fTrCorr;\r
+    iwd.ResetCovMat();\r
+    iwd.GetHitsPatt() = 0;\r
+    for (int i=pMax;i>=pLr;i--) {\r
+      KMCLayer* lr = GetLayer(i);\r
+      //      printf("IWD %d\n",i);\r
+      if (!lr->IsDead() && i!=pLr && i>=pMin) if (!UpdateTrack(&iwd,lr,&lr->fClCorr))  return 0;\r
+      if (i!=pLr) {\r
+       if (!iwd.CorrectForMeanMaterial(lr,kTRUE)) return 0; // correct for materials of this layer\r
+       if (!PropagateToLayer(&iwd,GetLayer(i-1),-1)) return 0;      // propagate to next layer\r
+      }\r
+      //  printf("IWD%d:  ",i); iwd.Print("l");\r
+    }\r
+    useIwd = kTRUE;\r
+  }\r
+  if (pLr>pMin) { // need outward piece\r
+    owd = GetLayer(pMin)->fTrCorr;\r
+    owd.ResetCovMat();\r
+    owd.GetHitsPatt() = 0;\r
+    for (int i=pMin;i<=pLr;i++) {\r
+      KMCLayer* lr = GetLayer(i);\r
+      //      printf("OWD %d\n",i);\r
+      if (!lr->IsDead() && i!=pLr && i<=pMax) if (!UpdateTrack(&owd,lr,&lr->fClCorr))  return 0;\r
+      if (i!=pLr) {\r
+       if (!owd.CorrectForMeanMaterial(lr,0)) return 0; // correct for materials of this layer\r
+       if (!PropagateToLayer(&owd,GetLayer(i+1), 1)) return 0;      // propagate to next layer\r
+      }\r
+      //      printf("OWD%d:  ",i); owd.Print("l");\r
+    }\r
+    useOwd = kTRUE;\r
+  }\r
+  //\r
+  // was this extrapolation outside the fit range?\r
+  if (!useIwd) return (KMCProbe*)&owd; \r
+  if (!useOwd) return (KMCProbe*)&iwd;\r
+  //\r
+  // weight both tracks\r
+  if (!iwd.Propagate(owd.GetAlpha(),owd.GetX(),fBFieldG)) return 0;\r
+  double meas[2] = {owd.GetY(),owd.GetZ()};\r
+  double measErr2[3] = {owd.GetSigmaY2(), owd.GetSigmaZY(), owd.GetSigmaZ2()};\r
+  //  printf("Weighting\n");\r
+  //  owd.Print("l");\r
+  //  iwd.Print("l");\r
+  if (!iwd.Update(meas,measErr2)) return 0;\r
+  iwd.GetHitsPatt() |= owd.GetHitsPatt();\r
+\r
+  //  printf("->\n");\r
+  //  iwd.Print("l");\r
+\r
+  return (KMCProbe*)&iwd;\r
+  //\r
+}\r
+\r
+//________________________________________________________________________________\r
+KMCProbe* KMCDetector::KalmanSmoothFull(int actLr, int actMin,int actMax) const\r
+{\r
+  // estimate kalman smoothed track params at given active lr \r
+  // from fit at layers actMin:actMax (excluding actLr)\r
+  // SolveSingleTrackViaKalman must have been called before\r
+  //\r
+  static TClonesArray prediction("KMCProbe",10);\r
+  static TClonesArray update("KMCProbe",10);\r
+  static KMCProbe res;\r
+  //\r
+  if (actMin>actMax) swap(actMin,actMax);\r
+  int nlrfit = actMax-actMin;\r
+  if (actLr>=actMin && actLr<=actMax) nlrfit-=1;\r
+  if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}\r
+  //\r
+  // find phisical layer id's\r
+  int pLr  = GetLayerID(actLr);\r
+  int pMin = GetLayerID(actMin);\r
+  int pMax = GetLayerID(actMax);\r
+  //\r
+  int dir=0,dirInt=0;\r
+  if      (pLr<=pMin) dir=-1; // inward extrapolation\r
+  else if (pLr>=pMax) dir= 1; // outward extrapolation\r
+  else if (actMax-actLr >= actLr-actMin) dirInt = -1; // inward  interpolation (the test point is closer to inner layer)\r
+  else    dirInt = 1;                                 // outward interpolation (the test point is closer to outer layer)\r
+  //\r
+  if (dir!=0) { // no sens to do smoothing: simple Kalman filtering extrapolation\r
+    int start = dir<0 ? pMax : pMin;\r
+    res = GetLayer(start)->fTrCorr;\r
+    res.ResetCovMat();\r
+    KMCLayer* lr = 0;\r
+    for (int i=(dir<0?pMax:pMin); i!=pLr; i+=dir) { // track till nearest layer to pLr\r
+      lr = GetLayer(i);\r
+      if (!lr->IsDead() && !(i<pMin ||i>pMax)) if (!UpdateTrack(&res,lr,&lr->fClCorr))  return 0; // update only with layers in fit range\r
+      if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this layer\r
+      if (!PropagateToLayer(&res,GetLayer(i+dir),dir))            return 0; // propagate to next layer     \r
+    }\r
+    if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this nearest layer\r
+    if (!PropagateToLayer(&res,GetLayer(pLr), dir)) return 0; // propagate to test layer\r
+    return (KMCProbe*)&res;\r
+  }\r
+  //\r
+  // too bad, need to do real filtering\r
+  //\r
+  int start = dirInt<0 ? pMax : pMin;\r
+  int stop  = dirInt<0 ? pMin-1 : pMax+1;\r
+  res = GetLayer(start)->fTrCorr;\r
+  res.ResetCovMat();\r
+  KMCLayer* lr = 0;\r
+  int count = 0;\r
+  for (int i=start; i!=stop; i+=dirInt) { // track in full range, storing updates and predictions\r
+    new(prediction[count]) KMCProbe(res);\r
+    lr = GetLayer(i);\r
+    if (!lr->IsDead() && i!=pLr) if (!UpdateTrack(&res,lr,&lr->fClCorr))  return 0; // update only with layers in fit range\r
+    new(update[count]) KMCProbe(res);\r
+    if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE))   return 0; // correct for materials of this layer\r
+    if (!PropagateToLayer(&res,GetLayer(i+dir),dir))            return 0; // propagate to next layer     \r
+    count++;\r
+  }\r
+  return (KMCProbe*)&res;\r
+  //\r
+}\r
+\r
+//________________________________________________________________________________\r
+Bool_t KMCDetector::SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta)\r
+{\r
+  // analytical estimate of tracking resolutions\r
+  //  fProbe.SetUseLogTermMS(kTRUE);\r
+  //\r
+  if (fMinITSHits>fNActiveITSLayers) {fMinITSHits = fNActiveITSLayers; printf("Redefined request of min N ITS hits to %d\n",fMinITSHits);}\r
+  if (TMath::Abs(eta)<1e-3) fDensFactorEta = 1.;\r
+  else {\r
+    fDensFactorEta = TMath::Tan( 2.*TMath::ATan(TMath::Exp(-TMath::Abs(eta))) );\r
+    fDensFactorEta = 1./TMath::Sqrt( 1. + 1./fDensFactorEta/fDensFactorEta);\r
+  }\r
+  double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta)); \r
+  KMCProbe* probe = PrepareKalmanTrack(pt,lambda,mass,-1);\r
+  if (!probe) return kFALSE;\r
+  //\r
+  KMCLayer *lr = 0;\r
+  //\r
+  //\r
+  // Start the track fitting --------------------------------------------------------\r
+  //\r
+  // Back-propagate the covariance matrix along the track. \r
+  // Kalman loop over the layers\r
+  //\r
+  KMCProbe* currTr = 0;\r
+  lr = (KMCLayer*)fLayers.At(fLastActiveLayerTracked);\r
+  lr->fTrCorr = *probe;\r
+  delete probe; // rethink...\r
+  //\r
+  for (Int_t j=fLastActiveLayerTracked; j--; ) {  // Layer loop\r
+    //\r
+    KMCLayer *lrP = lr;\r
+    lr = (KMCLayer*)fLayers.At(j);\r
+    //\r
+    lr->fTrCorr = lrP->fTrCorr;\r
+    currTr = &lr->fTrCorr;\r
+    currTr->ResetHit(lrP->GetActiveID());\r
+    //\r
+    // if there was a measurement on prev layer, update the track\r
+    if (!lrP->IsDead()) { // include measurement\r
+      KMCCluster cl(currTr->GetY(),currTr->GetZ(), currTr->GetX(), currTr->GetAlpha());\r
+      if (!UpdateTrack(currTr,lrP,&cl))  return kFALSE;\r
+    }\r
+    if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) return kFALSE; // correct for materials of this layer\r
+    if (!PropagateToLayer(currTr,lr,-1)) return kFALSE;      // propagate to current layer\r
+    //\r
+  } // end loop over layers\r
+  //\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________\r
+Bool_t KMCDetector::SolveSingleTrackViaKalmanMC(int offset)\r
+{\r
+  // MC estimate of tracking resolutions/effiencies. Requires that the SolveSingleTrackViaKalman\r
+  // was called before, since it uses data filled by this method\r
+  //\r
+  // The MC tracking will be done starting from fLastActiveITSLayer + offset (before analytical estimate will be used)\r
+  //\r
+  // At this point, the fProbe contains the track params generated at vertex.\r
+  // Clone it and propagate to target layer to generate hit positions affected by MS\r
+  //\r
+  fUpdCalls = 0.;\r
+  KMCProbe *currTrP=0,*currTr=0;\r
+  int maxLr = fLastActiveITSLayer + offset;\r
+  if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;\r
+  ResetMCTracks(maxLr);\r
+  KMCLayer* lr = (KMCLayer*)fLayers.At(maxLr);\r
+  currTr = lr->AddMCTrack(&fProbe); // start with original track at vertex\r
+  //\r
+  if (!TransportKalmanTrackWithMS(currTr, maxLr)) return kFALSE; // transport it to outermost layer where full MC is done\r
+  //\r
+  if (fLastActiveITSLayer<fLastActiveLayerTracked) { // prolongation from TPC\r
+    // start from correct track propagated from above till maxLr\r
+    double *covMS = (double*)currTr->GetCovariance();\r
+    const double *covIdeal =lr->fTrCorr.GetCovariance();\r
+    for (int i=15;i--;) covMS[i] = covIdeal[i];\r
+  }\r
+  else { // ITS SA: randomize the starting point\r
+    //    double *pars = (double*)currTr->GetParameter();\r
+    //    pars[0] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaY2()));\r
+    //    pars[1] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaZ2()));\r
+    //\r
+    currTr->ResetCovMat();\r
+    /*\r
+    double *trCov  = (double*)currTr->GetCovariance();\r
+    double *trPars = (double*)currTr->GetParameter();\r
+    const double kLargeErr2PtI = 0.3*0.3;\r
+    trCov[14] = TMath::Max(trCov[14],kLargeErr2PtI*trPars[4]*trPars[4]);\r
+    */\r
+  }\r
+  //\r
+  for (Int_t j=maxLr; j--; ) {  // Layer loop\r
+    //\r
+    KMCLayer *lrP = lr;\r
+    lr = (KMCLayer*)fLayers.At(j);\r
+    int ntPrev = lrP->GetNMCTracks();\r
+    //\r
+    if (lrP->IsDead()) { // for passive layer just propagate the copy of all tracks of prev layer >>>\r
+      for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
+       currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
+       currTr = lr->AddMCTrack( currTrP );\r
+       if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer\r
+       if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill(); continue;} // propagate to current layer\r
+      }\r
+      continue;\r
+    } // treatment of dead layer <<<\r
+    //\r
+    if (lrP->IsTPC()) { // we don't consider bg hits in TPC, just update with MC cluster\r
+      for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
+       currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
+       currTr = lr->AddMCTrack( currTrP );\r
+       if (!UpdateTrack(currTr, lrP, lrP->GetMCCluster(), kTRUE)) {currTr->Kill(); continue;} // update with correct MC cl.\r
+       if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer\r
+       if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill(); continue;} // propagate to current layer\r
+      }\r
+      continue;\r
+    } // treatment of ideal (TPC?) layer <<<\r
+    //\r
+    // active layer under eff. study (ITS?): propagate copy of every track to MC cluster frame (to have them all in the same frame)\r
+    // and calculate the limits of bg generation\r
+    KMCCluster* clMC = lrP->GetMCCluster();\r
+    if (lrP->GetLayerEff()<gRandom->Rndm()) clMC->Kill(); // simulate inefficiency\r
+    ResetSearchLimits();\r
+    int nseeds = 0;\r
+    for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
+      currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
+      currTr = lr->AddMCTrack( currTrP );\r
+      if (!currTr->PropagateToCluster(clMC,fBFieldG)) {currTr->Kill(); continue;} // go to MC cluster\r
+      if ( !(currTr->GetNITSHits()>0 && currTr->GetNITSHits()==currTr->GetNFakeITSHits()) ) UpdateSearchLimits(currTr, lrP); // RS\r
+      nseeds++;\r
+    }\r
+    //\r
+    //    printf("%3d seeds\n",nseeds);\r
+    if (fUseBackground && lrP->IsITS()) GenBgClusters(lrP); //  generate background hits\r
+    //\r
+    ntPrev = lr->GetNMCTracks();\r
+    for (int itr=ntPrev;itr--;) { // loop over all tracks PROPAGATED from previous layer to clusters frame on previous layer\r
+      currTrP = lr->GetMCTrack(itr); // this is a seed from prev layer. The new clusters are attached to its copies, the seed itself\r
+                                     // will be propagated w/o cluster update if it does not violate requested "reconstructed" track settings\r
+      if (currTrP->IsKilled()) continue;\r
+      //printf("Check    %d %p %d\n",itr,currTrP,currTrP->GetUniqueID()); currTrP->Print();\r
+      CheckTrackProlongations(currTrP, lr, lrP);\r
+      if (NeedToKill(currTrP)) currTrP->Kill(); // kill track which was not updated at lrP\r
+      //currTrP->Kill(); // kill track which was not updated at lrP\r
+    }\r
+    //  \r
+    lr->GetMCTracks()->Sort();\r
+    int ntTot = lr->GetNMCTracks(); // propagate max amount of allowed tracks to current layer\r
+    if (ntTot>fMaxSeedToPropagate && fMaxSeedToPropagate>0) {\r
+      for (int itr=ntTot;itr>=fMaxSeedToPropagate;itr--)  lr->GetMCTracks()->RemoveAt(itr);\r
+      ntTot = fMaxSeedToPropagate;\r
+    }\r
+    //\r
+    for (int itr=ntTot;itr--;) {\r
+      currTr = lr->GetMCTrack(itr);\r
+      if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill();continue;} // correct for materials of prev. layer\r
+      if (!PropagateToLayer(currTr,lr,-1))      {currTr->Kill();continue;} // propagate to current layer\r
+    }\r
+    AliDebug(1,Form("Got %d tracks on layer %s",ntTot,lr->GetName()));\r
+    //    lr->GetMCTracks()->Print();\r
+    //\r
+  } // end loop over layers    \r
+  //\r
+  // do we use vertex constraint?\r
+  KMCLayer *vtx = GetLayer(0);\r
+  if (!vtx->IsDead() && vtx->IsITS()) {\r
+    int ntr = vtx->GetNMCTracks();\r
+    for (int itr=0;itr<ntr;itr++) {\r
+      currTr = vtx->GetMCTrack(itr);\r
+      if (currTr->IsKilled()) continue;\r
+      KMCCluster* clv = vtx->GetMCCluster();\r
+      double meas[2] = {clv->GetY(),clv->GetZ()};\r
+      double measErr2[3] = {vtx->fPhiRes*vtx->fPhiRes,0,vtx->fZRes*vtx->fZRes};\r
+      double chi2v = currTr->GetPredictedChi2(meas,measErr2);\r
+      currTr->AddHit(vtx->GetActiveID(), chi2v, -1);\r
+      currTr->SetInnerLrChecked(vtx->GetActiveID());\r
+      if (NeedToKill(currTr)) currTr->Kill();\r
+      // if (vtx->IsITS()) {if (!UpdateTrack(currTr, vtx, vtx->GetMCCluster(), kFALSE)) {currTr->Kill();continue;}}\r
+    }\r
+  }\r
+  EliminateUnrelated();\r
+  \r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t KMCDetector::PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const\r
+{\r
+  // bring the track to layer and rotat to frame normal to its surface\r
+  if (!trc->PropagateToR(lr->fR,fBFieldG, dir)) return kFALSE;\r
+  //\r
+  // rotate to frame with X axis normal to the surface (defined by ideal track)\r
+  if (!lr->IsVertex()) {\r
+    double pos[3];\r
+    trc->GetXYZ(pos);  // lab position\r
+    double phi = TMath::ATan2(pos[1],pos[0]);\r
+    if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;\r
+    if (!trc->Rotate(phi)) {\r
+      AliDebug(2,Form("Failed to rotate to the frame (phi:%+.3f) of layer at %.2f at XYZ: %+.3f %+.3f %+.3f",\r
+                     phi,lr->fR,pos[0],pos[1],pos[2]));\r
+      if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");\r
+      return kFALSE;\r
+    }\r
+  }\r
+  //\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t KMCDetector::UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster) const\r
+{\r
+  // update track with measured cluster\r
+  // propagate to cluster\r
+  double meas[2] = {cl->GetY(),cl->GetZ()}; // ideal cluster coordinate\r
+  double measErr2[3] = {lr->fPhiRes*lr->fPhiRes,0,lr->fZRes*lr->fZRes};\r
+  //\r
+  if (goToCluster) if (!trc->PropagateToCluster(cl,fBFieldG)) return kFALSE; // track was not propagated to cluster frame\r
+  //\r
+  double chi2 = trc->GetPredictedChi2(meas,measErr2);\r
+  //  if (chi2>fMaxChi2Cl) return kTRUE; // chi2 is too large\r
+  //  \r
+  if (!trc->Update(meas,measErr2)) {\r
+    AliDebug(2,Form("layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",\r
+                   lr->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));\r
+    if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");\r
+    return kFALSE;\r
+  }\r
+  trc->AddHit(lr->GetActiveID(), chi2);\r
+  //\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________________________\r
+Int_t KMCDetector::GenBgClusters(KMCLayer* lr)\r
+{\r
+  // Generate fake clusters in precalculated RPhi,Z range\r
+  if (fNBgLimits<1) return 0; // limits were not set - no track was prolongated\r
+  //\r
+  // Fix search limits to avoid seeds which will anyway point very far from the vertex\r
+  double tolY = TMath::Sqrt(lr->fSig2EstD)*fMaxChi2ClSQ;\r
+  double tolZ = TMath::Sqrt(lr->fSig2EstZ)*fMaxChi2ClSQ;\r
+  \r
+  //  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);\r
+  if (fBgYMin < lr->fClCorr.fY-tolY) fBgYMin = lr->fClCorr.fY-tolY;\r
+  if (fBgYMax > lr->fClCorr.fY+tolY) fBgYMax = lr->fClCorr.fY+tolY;\r
+  if (fBgZMin < lr->fClCorr.fZ-tolZ) fBgZMin = lr->fClCorr.fZ-tolZ;\r
+  if (fBgZMax > lr->fClCorr.fZ+tolZ) fBgZMax = lr->fClCorr.fZ+tolZ;\r
+  //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);\r
+  //\r
+  double dy = fBgYMax - fBgYMin;\r
+  double dz = fBgZMax - fBgZMin;\r
+  double surf = dy*dz;               // surface of generation\r
+  if (surf<0) return 0;\r
+  double poissProb = surf*HitDensity(lr->fR)*lr->GetLayerEff();\r
+  AliDebug(2,Form("Bg for Lr %s (r=%.2f) : Density %.2f on surface %.2e [%+.4f : %+.4f][%+.4f %+.4f]",\r
+                 lr->GetName(),lr->fR,HitDensity(lr->fR),surf,fBgYMin,fBgYMax,fBgZMin,fBgZMax));\r
+  int nFakesGen = gRandom->Poisson( poissProb ); // preliminary number of extra clusters to test\r
+  KMCCluster *refCl = lr->GetMCCluster();\r
+  double sig2y = lr->GetPhiRes()*lr->GetPhiRes();\r
+  double sig2z = lr->GetZRes()*lr->GetZRes();\r
+  for (int ic=nFakesGen;ic--;) {\r
+    double y = fBgYMin+dy*gRandom->Rndm();\r
+    double z = fBgZMin+dz*gRandom->Rndm();\r
+    double dfy = y-refCl->GetY();\r
+    double dfz = z-refCl->GetZ();\r
+    double dist = (dfy*dfy)/sig2y + (dfz*dfz)/sig2z;\r
+    if (dist<4) continue; // avoid overlap with MC cluster\r
+    lr->AddBgCluster(y, z, refCl->GetX(), refCl->GetPhi());\r
+  }\r
+  AliDebug(2,Form("Added %6d noise clusters on lr %s (poisson Prob=%8.2f for surface %.2e) DY:%7.4f DZ: %7.4f",\r
+                  lr->GetNBgClusters(),lr->GetName(),poissProb,surf,dy,dz));\r
+  return nFakesGen;\r
+  //\r
+}\r
+\r
+//____________________________________________________________________________\r
+void KMCDetector::UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr)\r
+{\r
+  // define the search window for track on layer (where the bg hist will be generated)\r
+  static double *currYMin = fBgYMinTr.GetArray();\r
+  static double *currYMax = fBgYMaxTr.GetArray();\r
+  static double *currZMin = fBgZMinTr.GetArray();\r
+  static double *currZMax = fBgZMaxTr.GetArray();\r
+  //\r
+  double sizeY = probe->GetSigmaY2(), sizeZ = probe->GetSigmaZ2();\r
+  //  /*\r
+  if (probe->GetNITSHits()<3) sizeY = 10*lr->fSig2EstD;\r
+  if (probe->GetNITSHits()<2) sizeZ = 10*lr->fSig2EstZ;\r
+  sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY+lr->fPhiRes*lr->fPhiRes); // max deviation in rphi to accept\r
+  sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ+lr->fZRes*lr->fZRes); // max deviation in dz to accept\r
+  //  */\r
+  //\r
+  /*\r
+  if (probe->GetNITSHits()<3) sizeY = 1./(1./lr->fSig2EstD + 1./sizeY);\r
+  if (probe->GetNITSHits()<2) sizeZ = 1./(1./lr->fSig2EstZ + 1./sizeZ);\r
+  sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY); // max deviation in rphi to accept\r
+  sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ); // max deviation in dz to accept\r
+  */\r
+  //\r
+  //  if (sizeY>2) sizeY=2;\r
+  //  if (sizeZ>2) sizeZ=2;\r
+  //  printf("Sizes at %s: %.5f %.5f\n",lr->GetName(), sizeY,sizeZ);\r
+  //\r
+  if (fNBgLimits>=fBgYMinTr.GetSize()) { // expand arrays, update pointers\r
+    fBgYMinTr.Set(2*(fNBgLimits+1));\r
+    fBgYMaxTr.Set(2*(fNBgLimits+1));\r
+    fBgZMinTr.Set(2*(fNBgLimits+1));\r
+    fBgZMaxTr.Set(2*(fNBgLimits+1));\r
+    currYMin = fBgYMinTr.GetArray();\r
+    currYMax = fBgYMaxTr.GetArray();\r
+    currZMin = fBgZMinTr.GetArray();\r
+    currZMax = fBgZMaxTr.GetArray();\r
+  }\r
+  if (fBgYMin > (currYMin[fNBgLimits]=probe->GetY()-sizeY) ) fBgYMin = currYMin[fNBgLimits];\r
+  if (fBgYMax < (currYMax[fNBgLimits]=probe->GetY()+sizeY) ) fBgYMax = currYMax[fNBgLimits];\r
+  if (fBgZMin > (currZMin[fNBgLimits]=probe->GetZ()-sizeZ) ) fBgZMin = currZMin[fNBgLimits];\r
+  if (fBgZMax < (currZMax[fNBgLimits]=probe->GetZ()+sizeZ) ) fBgZMax = currZMax[fNBgLimits];\r
+  if (AliLog::GetGlobalDebugLevel()>=2) {\r
+    probe->Print("l");\r
+    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]));\r
+    AliInfo(Form("Global Limits Lr %s                            [%+.4f : %+.4f][%+.4f %+.4f]",lr->GetName(),fBgYMin,fBgYMax,fBgZMin,fBgZMax));\r
+    AliInfo(Form("MC Cluster: %+.4f : %+.4f",lr->fClMC.fY, lr->fClMC.fZ));\r
+  }\r
+  probe->SetUniqueID(fNBgLimits++);\r
+  //\r
+  if (lr->IsITS() && probe->GetNFakeITSHits()==0) {\r
+    if (fHMCLrResidRPhi) fHMCLrResidRPhi->Fill(probe->GetY() - lr->GetMCCluster()->GetY(), lr->GetActiveID());\r
+    if (fHMCLrResidZ)    fHMCLrResidZ->Fill(probe->GetZ() - lr->GetMCCluster()->GetZ(),lr->GetActiveID());\r
+  }\r
+  //\r
+}\r
+\r
+//____________________________________________________________________________\r
+void KMCDetector::CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP)\r
+{\r
+  // explore prolongation of probe from lrP to lr with all possible clusters of lrP\r
+  // the probe is already brought to clusters frame\r
+  int nCl = lrP->GetNBgClusters();\r
+  double measErr2[3] = {lrP->fPhiRes*lrP->fPhiRes,0,lrP->fZRes*lrP->fZRes};\r
+  double meas[2] = {0,0};\r
+  UInt_t tmpID = probe->GetUniqueID();\r
+  double yMin = fBgYMinTr[tmpID];\r
+  double yMax = fBgYMaxTr[tmpID];\r
+  double zMin = fBgZMinTr[tmpID];\r
+  double zMax = fBgZMaxTr[tmpID];\r
+  //\r
+  probe->SetInnerLrChecked(lrP->GetActiveID());\r
+  for (int icl=-1;icl<nCl;icl++) {\r
+    KMCCluster* cl = icl<0 ? lrP->GetMCCluster() : lrP->GetBgCluster(icl);  // -1 is for true MC cluster\r
+    if (cl->IsKilled()) {\r
+      if (AliLog::GetGlobalDebugLevel()>1) {printf("Skip cluster %d ",icl); cl->Print();}\r
+      continue;\r
+    }\r
+    double y = cl->GetY();\r
+    double z = cl->GetZ();\r
+    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));\r
+    if (AliLog::GetGlobalDebugLevel()>0) {\r
+      if (icl==-1 && probe->GetNFakeITSHits()==0) {\r
+       meas[0] = y; meas[1] = z;\r
+       double chi2a = probe->GetPredictedChi2(meas,measErr2);\r
+       if (chi2a>fMaxChi2Cl || (y<yMin || y>yMax) || (z<zMin || z>zMax)) {\r
+         probe->Print();\r
+         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",\r
+                y,z,lrP->GetName(),chi2a,y-probe->GetY(),z-probe->GetZ(),yMin,yMax,zMin,zMax, probe->GetX(), cl->GetX(), probe->GetAlpha(), cl->GetPhi());\r
+       }\r
+      }\r
+    }\r
+    if (y<yMin || y>yMax) continue; // preliminary check on Y\r
+    if (z<zMin || z>zMax) continue; // preliminary check on Z\r
+    meas[0] = y; meas[1] = z;\r
+    double chi2 = probe->GetPredictedChi2(meas,measErr2);\r
+    if (fHMCLrChi2 && probe->GetNFakeITSHits()==0 && icl==-1) fHMCLrChi2->Fill(chi2,lrP->GetActiveID());\r
+    AliDebug(2,Form("Seed-to-cluster chi2 = Chi2=%.2f",chi2));\r
+    if (chi2>fMaxChi2Cl) continue;\r
+    // \r
+    // update track copy\r
+    KMCProbe* newTr = lr->AddMCTrack( probe );\r
+    fUpdCalls++;\r
+    if (!newTr->Update(meas,measErr2)) {\r
+      AliDebug(2,Form("Layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",\r
+                     lrP->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));\r
+      if (AliLog::GetGlobalDebugLevel()>1) newTr->Print("l");\r
+      newTr->Kill();\r
+      continue;\r
+    }\r
+    newTr->AddHit(lrP->GetActiveID(), chi2, icl);\r
+    if (AliLog::GetGlobalDebugLevel()>1) {\r
+      AliInfo("Cloned updated track is:");\r
+      newTr->Print();\r
+    }\r
+    if (NeedToKill(newTr)) newTr->Kill();\r
+  }\r
+  //\r
+}\r
+\r
+//____________________________________________________________________________\r
+void KMCDetector::ResetMCTracks(Int_t maxLr)\r
+{\r
+  int nl = GetNLayers();\r
+  if (maxLr<0 || maxLr>=nl) maxLr = nl-1;\r
+  for (int i=maxLr+1;i--;) GetLayer(i)->ResetMCTracks();\r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t KMCDetector::NeedToKill(KMCProbe* probe) const\r
+{\r
+  // check if the seed at given layer (last one where update was tried) \r
+  // still has chances to be reconstructed\r
+  const Bool_t kModeKillMiss = kFALSE;\r
+  //\r
+  Bool_t kill = kFALSE;\r
+  while (1) {\r
+    int il = probe->GetInnerLayerChecked();\r
+    int nITS = probe->GetNITSHits();\r
+    int nITSMax = nITS + il; // maximum it can have\r
+    if (nITSMax<fMinITSHits) {kill = kTRUE; break;}    // has no chance to collect enough ITS hits\r
+    //\r
+    int ngr = fPattITS.GetSize();\r
+    if (ngr>0) { // check pattern\r
+      UInt_t patt = probe->GetHitsPatt();\r
+      // complete the layers not checked yet\r
+      for (int i=il;i--;) patt |= (0x1<<i);\r
+      for (int ig=ngr;ig--;) if (!(((UInt_t)fPattITS[ig]) & patt)) {kill = kTRUE; break;}\r
+      //\r
+    }\r
+    //\r
+    if (nITS>2) {  // check if smallest possible norm chi2/ndf is acceptable\r
+      double chi2min = probe->GetChi2();\r
+      if (kModeKillMiss) {\r
+       int nMiss = fNActiveITSLayers - probe->GetInnerLayerChecked() - nITS; // layers already missed\r
+       chi2min = nMiss*probe->GetMissingHitPenalty();\r
+      }\r
+      chi2min /= ((nITSMax<<1)-KMCProbe::kNDOF);\r
+      if (chi2min>fMaxNormChi2NDF) {kill = kTRUE; break;}\r
+    }\r
+    //\r
+    // loose vertex constraint\r
+    double dst;\r
+    if (nITS>=2) {\r
+      probe->GetZAt(0,fBFieldG,dst);\r
+      //printf("Zd (F%d): %f\n",probe->GetNFakeITSHits(),dst);\r
+      if (TMath::Abs(dst)>10.) {kill = kTRUE; break;}\r
+    }\r
+    if (nITS>=3) {\r
+      probe->GetYAt(0,fBFieldG,dst);\r
+      //printf("Dd (F%d): %f\n",probe->GetNFakeITSHits(),dst);\r
+      if (TMath::Abs(dst)>10.) {kill = kTRUE; break;}\r
+    }\r
+    //\r
+    break;\r
+  }\r
+  if (kill && AliLog::GetGlobalDebugLevel()>1 && probe->GetNFakeITSHits()==0) {\r
+    printf("Killing good seed, last upd layer was %d\n",probe->GetInnerLayerChecked());\r
+    probe->Print("l");\r
+  }\r
+  return kill;\r
+}\r
+\r
+//_____________________________________________________________________\r
+void KMCDetector::EliminateUnrelated()\r
+{\r
+  // kill useless tracks\r
+  KMCLayer* lr = GetLayer(0);\r
+  int ntr = lr->GetNMCTracks();\r
+  int nval = 0;\r
+  for (int itr=0;itr<ntr;itr++) {\r
+    KMCProbe* probe = lr->GetMCTrack(itr);\r
+    if (probe->IsKilled()) continue;\r
+    if (probe->GetNITSHits()-probe->GetNFakeITSHits()<1) {probe->Kill(); continue;}\r
+    nval++;\r
+  }\r
+  lr->GetMCTracks()->Sort();\r
+  const int kDump = 0;\r
+  if (kDump>0) {\r
+    printf("Valid %d out of %d\n",nval, ntr);\r
+    ntr = ntr>kDump ? kDump:0;\r
+    for (int itr=0;itr<ntr;itr++) {\r
+      lr->GetMCTrack(itr)->Print();\r
+    }\r
+  }\r
+}\r
+\r
+//_____________________________________________________________________\r
+void KMCDetector::RequirePattern(UInt_t *patt, int groups)\r
+{\r
+  if (groups<1) {fPattITS.Set(0); return;}\r
+  fPattITS.Set(groups);\r
+  for (int i=0;i<groups;i++) fPattITS[i] = patt[i];\r
+}\r
+\r
+//_____________________________________________________________________\r
+void KMCDetector::CalcHardSearchLimits(double dzv)\r
+{\r
+  //\r
+  TArrayD zlims;\r
+  zlims.Set(fNActiveITSLayers);\r
+  for (int il0=0;il0<fNActiveITSLayers;il0++) {\r
+    KMCLayer* lr0 = GetActiveLayer(il0);\r
+    double angZTol = dzv/lr0->GetRadius();\r
+    for (int il1=0;il1<fNActiveITSLayers;il1++) {\r
+      if (il1==il0) continue;\r
+      KMCLayer* lr1 = GetActiveLayer(il1);\r
+      double ztol = angZTol*TMath::Abs(lr0->GetRadius() - lr1->GetRadius());\r
+      if (ztol>zlims[il1]) zlims[il1] = ztol;\r
+    }\r
+  }\r
+  //\r
+  for (int il=0;il<fNActiveITSLayers;il++) printf("ZTol%d: %8.4f\n",il,zlims[il]);\r
+}\r
+\r
+//_______________________________________________________\r
+double KMCDetector::PropagateBack(KMCProbe* trc) \r
+{\r
+  static KMCProbe bwd;\r
+  bwd = *trc;\r
+  bwd.ResetCovMat();\r
+  static double measErr2[3] = {0,0,0};\r
+  static double meas[2] = {0,0};\r
+  int icl = 0;\r
+  double chi2Tot = 0;\r
+  for (int il=1;il<=fLastActiveITSLayer;il++) {\r
+    KMCLayer* lr = GetLayer(il);\r
+    if (!PropagateToLayer(&bwd,lr,1)) return -1;\r
+    int aID = lr->GetActiveID();\r
+    if (aID>-1 && (icl=bwd.fClID[aID])>=-1) {\r
+      KMCCluster* clMC =  icl<0 ? lr->GetMCCluster() : lr->GetBgCluster(icl);\r
+      if (!bwd.PropagateToCluster(clMC,fBFieldG)) return -1;\r
+      meas[0] = clMC->GetY(); meas[1] = clMC->GetZ();\r
+      measErr2[0] = lr->fPhiRes*lr->fPhiRes;\r
+      measErr2[2] = lr->fZRes*lr->fZRes;\r
+      double chi2a = bwd.GetPredictedChi2(meas,measErr2);\r
+      chi2Tot += chi2a;\r
+      printf("Chis %d (cl%+3d): t2c: %6.3f tot: %6.3f\n",aID,icl,chi2a, chi2Tot);\r
+      bwd.Update(meas,measErr2);\r
+      bwd.AddHit(aID, chi2a, icl);     \r
+    }\r
+    if (!bwd.CorrectForMeanMaterial(lr,kFALSE)) return -1;\r
+  }\r
+  return chi2Tot;\r
+}\r
+\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+ClassImp(KMCTrackSummary)\r
+\r
+Int_t KMCTrackSummary::fgSumCounter = 0;\r
+\r
+//____________________________________________________________________________\r
+KMCTrackSummary::KMCTrackSummary(const char* name,const char* title, int nlr) :\r
+  TNamed(name,name), fNITSLayers(nlr),\r
+  fPattITS(0),fPattITSFakeExcl(0),fPattITSCorr(0),\r
+  fHMCChi2(0),fHMCSigDCARPhi(0),fHMCSigDCAZ(0),fHMCSigPt(0),fHMCNCl(0),fHMCFakePatt(0),\r
+  fCountAcc(0),fCountTot(0),fUpdCalls(0),\r
+  fRefProbe(),fAnProbe()\r
+{\r
+  // create summary structure for specific track conditions\r
+  //\r
+  SetMinMaxClITS();\r
+  SetMinMaxClITSFake();\r
+  SetMinMaxClITSCorr();\r
+  //\r
+  if (name==0 && title==0 && nlr==0) return;  // default dummy contructor\r
+  //\r
+  enum {kNBinRes=1000};\r
+  const double kMaxChi2(10),kMaxResRPhi=0.1, kMaxResZ=0.1,kMaxResPt=0.3;\r
+  //\r
+  TString strN = name, strT = title;\r
+  if (strN.IsNull()) strN = Form("TrSum%d",fgSumCounter);\r
+  if (strT.IsNull()) strT = strN;\r
+  fgSumCounter++;\r
+  //\r
+  if (fNITSLayers<1) {\r
+    fNITSLayers=KMCProbe::GetNITSLayers(); \r
+    if (fNITSLayers<1) {AliError("N ITS layer is not provided and not available from KMCProbe::GetNITSLayers()"); exit(1);}\r
+    AliInfo(Form("%s No N Layers provided, set %d",strN.Data(),fNITSLayers));\r
+  }\r
+  nlr = fNITSLayers;\r
+  //\r
+  TString nam,tit;\r
+  //\r
+  nam = Form("%s_mc_chi2",strN.Data());\r
+  tit = Form("%s MC #chi^{2}",strT.Data());\r
+  fHMCChi2 = new TH1F(nam.Data(),tit.Data(),kNBinRes,0,kMaxChi2);\r
+  fHMCChi2->GetXaxis()->SetTitle("#chi^{2}");\r
+  fHMCChi2->Sumw2();\r
+  //\r
+  nam = Form("%s_mc_DCArphi",strN.Data());\r
+  tit = Form("%s MC #DeltaDCA r#phi",strT.Data());\r
+  fHMCSigDCARPhi = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResRPhi,kMaxResRPhi);\r
+  fHMCSigDCARPhi->GetXaxis()->SetTitle("#Delta r#phi");\r
+  fHMCSigDCARPhi->Sumw2();\r
+  //\r
+  nam = Form("%s_mc_DCAz",strN.Data());\r
+  tit = Form("%s MC #DeltaDCA Z",strT.Data());\r
+  fHMCSigDCAZ = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResZ,kMaxResZ);\r
+  fHMCSigDCAZ->GetXaxis()->SetTitle("#Delta Z");\r
+  fHMCSigDCAZ->Sumw2();\r
+  //\r
+  nam = Form("%s_mc_pt",strN.Data());\r
+  tit = Form("%s MC $Deltap_{T}/p_{T}",strT.Data());\r
+  fHMCSigPt = new TH1F(nam.Data(),tit.Data(),2*kNBinRes,-kMaxResPt,kMaxResPt);\r
+  fHMCSigPt->GetXaxis()->SetTitle("#Deltap_{T}/p_{T}");\r
+  fHMCSigPt->Sumw2();\r
+  //\r
+  nam = Form("%s_n_cl",strN.Data());\r
+  tit = Form("%s N Clusters",strT.Data());\r
+  fHMCNCl = new TH2F(nam.Data(),tit.Data(),nlr,1,nlr+1,nlr,1,nlr+1);\r
+  fHMCNCl->GetXaxis()->SetTitle("N Clusters Tot");\r
+  fHMCNCl->GetYaxis()->SetTitle("N Clusters Fake");\r
+  fHMCNCl->Sumw2();\r
+  //\r
+  nam = Form("%s_fake_cl",strN.Data());\r
+  tit = Form("%s Fake Clusters",strT.Data());\r
+  fHMCFakePatt = new TH1F(nam.Data(),tit.Data(),nlr,-0.5,nlr-0.5);\r
+  fHMCFakePatt->GetXaxis()->SetTitle("Fake clusters pattern");\r
+  fHMCFakePatt->Sumw2();\r
+  //\r
+}\r
+\r
+//____________________________________________________________________________\r
+KMCTrackSummary::~KMCTrackSummary()\r
+{\r
+  if (fHMCChi2)         delete fHMCChi2;\r
+  if (fHMCSigDCARPhi)   delete fHMCSigDCARPhi;\r
+  if (fHMCSigDCAZ)      delete fHMCSigDCAZ;\r
+  if (fHMCSigPt)        delete fHMCSigPt;\r
+  if (fHMCNCl)          delete fHMCNCl;\r
+  if (fHMCFakePatt)     delete fHMCFakePatt;\r
+}\r
+\r
+\r
+//____________________________________________________________________________\r
+void KMCTrackSummary::Add(const KMCTrackSummary* src, double scl)\r
+{\r
+  if (fHMCChi2)         fHMCChi2->Add(src->fHMCChi2,scl);\r
+  if (fHMCSigDCARPhi)   fHMCSigDCARPhi->Add(src->fHMCSigDCARPhi,scl);\r
+  if (fHMCSigDCAZ)      fHMCSigDCAZ->Add(src->fHMCSigDCAZ,scl);\r
+  if (fHMCSigPt)        fHMCSigPt->Add(src->fHMCSigPt,scl);\r
+  if (fHMCNCl)          fHMCNCl->Add(src->fHMCNCl,scl);\r
+  if (fHMCFakePatt)     fHMCFakePatt->Add(src->fHMCFakePatt,scl);\r
+  fCountAcc += src->fCountAcc;\r
+  fUpdCalls += src->fUpdCalls;\r
+}\r
+\r
+//____________________________________________________________________________\r
+void KMCTrackSummary::PrependTNamed(TNamed* obj, const char* nm, const char* tit)\r
+{\r
+  // prepend  name, title by this prefix\r
+  TString name = obj->GetName(),title = obj->GetTitle();\r
+  name.Prepend(nm); title.Prepend(tit);\r
+  obj->SetNameTitle(name.Data(),title.Data());\r
+  //\r
+}\r
+  \r
+//____________________________________________________________________________\r
+void KMCTrackSummary::SetNamePrefix(const char* pref)\r
+{\r
+  // prepend all names by this prefix\r
+  TString nmT,nmP = pref;\r
+  if (nmP.IsNull()) return;\r
+  nmP = Form("%s_",pref);\r
+  nmT = Form("%s ",pref);\r
+  PrependTNamed(this, nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCChi2,       nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCSigDCARPhi, nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCSigDCAZ,    nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCSigPt,      nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCNCl,        nmP.Data(), nmT.Data());\r
+  PrependTNamed(fHMCFakePatt,   nmP.Data(), nmT.Data());\r
+  //  \r
+}\r
+\r
+//____________________________________________________________________________\r
+Bool_t KMCTrackSummary::CheckTrack(KMCProbe* trc)\r
+{\r
+  // check if the track satisfies to selection conditions\r
+  fCountTot++;\r
+  if (!trc) return kFALSE;\r
+  if (OutOfRange(trc->GetNITSHits(), fMinMaxClITS)) return kFALSE;\r
+  if (OutOfRange(trc->GetNFakeITSHits(), fMinMaxClITSFake)) return kFALSE;\r
+  if (OutOfRange(trc->GetNITSHits()-trc->GetNFakeITSHits(), fMinMaxClITSCorr)) return kFALSE;\r
+  //\r
+  // check layer patterns if requested\r
+  UInt_t patt  = trc->GetHitsPatt();\r
+  UInt_t pattF = trc->GetFakesPatt();\r
+  UInt_t pattC = patt&(~pattF);    // correct hits\r
+  // is there at least one hit in each of requested layer groups?\r
+  for (int ip=fPattITS.GetSize();ip--;)     if (!CheckPattern(patt,fPattITS[ip])) return kFALSE;\r
+  // is there at least one fake in any of requested layer groups?\r
+  for (int ip=fPattITSFakeExcl.GetSize();ip--;) if (CheckPattern(pattF,fPattITSFakeExcl[ip])) return kFALSE;\r
+  // is there at least one correct hit in each of requested layer groups?\r
+  for (int ip=fPattITSCorr.GetSize();ip--;) if (!CheckPattern(pattC,fPattITSCorr[ip])) return kFALSE;\r
+  //\r
+  fCountAcc++;\r
+  return kTRUE;\r
+}\r
+\r
+//____________________________________________________________________________\r
+void KMCTrackSummary::AddTrack(KMCProbe* trc)\r
+{\r
+  // fill track info\r
+  if (!CheckTrack(trc)) return;\r
+  //\r
+  fHMCChi2->Fill(trc->GetNormChi2(kFALSE));\r
+  fHMCSigDCARPhi->Fill(trc->GetY());\r
+  fHMCSigDCAZ->Fill(trc->GetZ());\r
+  if (fRefProbe.Pt()>0) fHMCSigPt->Fill( trc->Pt()/fRefProbe.Pt()-1.);\r
+  //  printf("Pts: %.3f %.3f -> %.3f\n",trc->Pt(),fRefProbe.Pt(),trc->Pt()/fRefProbe.Pt()-1.);\r
+  //  trc->Print("l");\r
+  //  fRefProbe.Print("l");\r
+  fHMCNCl->Fill(trc->GetNITSHits(),trc->GetNFakeITSHits());\r
+  for (int i=trc->GetNITSLayers();i--;) if (trc->IsHitFake(i)) fHMCFakePatt->Fill(i);\r
+  //\r
+}\r
+\r
+//____________________________________________________________________________\r
+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,\r
+                         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,\r
+                         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,\r
+                         Bool_t l30,Bool_t l31)\r
+{\r
+  // create corresponding bit pattern\r
+  UInt_t patt = 0;\r
+  if (l0 ) patt |= (0x1<<0);    if (l1 ) patt |= (0x1<<1);   if (l2 ) patt |= (0x1<<2);\r
+  if (l3 ) patt |= (0x1<<3);    if (l4 ) patt |= (0x1<<4);   if (l5 ) patt |= (0x1<<5);\r
+  if (l6 ) patt |= (0x1<<6);    if (l7 ) patt |= (0x1<<7);   if (l8 ) patt |= (0x1<<8);\r
+  if (l9 ) patt |= (0x1<<9);    if (l10) patt |= (0x1<<10);  if (l11) patt |= (0x1<<11);\r
+  if (l12) patt |= (0x1<<12);   if (l13) patt |= (0x1<<13);  if (l14) patt |= (0x1<<14);\r
+  if (l15) patt |= (0x1<<15);   if (l16) patt |= (0x1<<16);  if (l17) patt |= (0x1<<17);\r
+  if (l18) patt |= (0x1<<18);   if (l19) patt |= (0x1<<19);  if (l20) patt |= (0x1<<20);\r
+  if (l21) patt |= (0x1<<21);   if (l22) patt |= (0x1<<22);  if (l23) patt |= (0x1<<23);\r
+  if (l24) patt |= (0x1<<24);   if (l25) patt |= (0x1<<25);  if (l26) patt |= (0x1<<26);\r
+  if (l27) patt |= (0x1<<27);   if (l28) patt |= (0x1<<28);  if (l29) patt |= (0x1<<29);\r
+  if (l30) patt |= (0x1<<30);   if (l31) patt |= (0x1<<31);   \r
+  return patt;\r
+}\r
+\r
+//__________________________________________________________________________\r
+void KMCTrackSummary::Print(Option_t* ) const\r
+{\r
+  printf("%s: summary for track M=%5.3f pT: %6.3f eta: %.2f\n", GetName(),\r
+        fRefProbe.GetMass(),fRefProbe.Pt(), fRefProbe.Eta());\r
+  printf("Cuts on NCl ITS: Tot: %2d - %2d Fake: %2d - %2d Corr: %2d - %2d\n",\r
+        fMinMaxClITS[0],fMinMaxClITS[1], \r
+        fMinMaxClITSFake[0],fMinMaxClITSFake[1],\r
+        fMinMaxClITSCorr[0],fMinMaxClITSCorr[1]);\r
+  //\r
+  int nlr = fNITSLayers;\r
+  if (fPattITS.GetSize()) {\r
+    printf("Require at least 1 hit in groups: ");\r
+    printf("Hits obligatory in groups: ");\r
+    for (int i=fPattITS.GetSize();i--;) {\r
+      UInt_t pat = (UInt_t)fPattITS[i];\r
+      printf("[");\r
+      for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
+      printf("] ");\r
+    }\r
+    printf("\n");\r
+  }\r
+  //\r
+  if (fPattITSFakeExcl.GetSize()) {\r
+    printf("Fake hit forbidden in groups    : ");\r
+    for (int i=fPattITSFakeExcl.GetSize();i--;) {\r
+      UInt_t pat = (UInt_t)fPattITSFakeExcl[i];\r
+      printf("[");\r
+      for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
+      printf("] ");\r
+    }\r
+    printf("\n");\r
+  }\r
+  //\r
+  if (fPattITSCorr.GetSize()) {\r
+    printf("Correct hit obligatory in groups: ");\r
+    for (int i=fPattITSCorr.GetSize();i--;) {\r
+      UInt_t pat = (UInt_t)fPattITSCorr[i];\r
+      printf("[");\r
+      for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
+      printf("] ");\r
+    }\r
+    printf("\n");\r
+  }\r
+  //  \r
+  printf("Entries: Tot: %.4e Acc: %.4e -> Eff: %.4f+-%.4f %.2e KMC updates\n",fCountTot,fCountAcc,GetEff(),GetEffErr(),GetUpdCalls());\r
+}\r
+\r
diff --git a/ITS/UPGRADE/KMCDetector.h b/ITS/UPGRADE/KMCDetector.h
new file mode 100755 (executable)
index 0000000..f66b528
--- /dev/null
@@ -0,0 +1,550 @@
+#ifndef DETECTORK_H\r
+#define DETECTORK_H\r
+\r
+#include <TNamed.h>\r
+#include <TList.h>\r
+#include <TGraph.h>\r
+#include <Riostream.h>\r
+#include <TMatrixD.h>\r
+#include <TClonesArray.h>\r
+#include <TArrayD.h>\r
+#include <TH2.h>\r
+#include <TH1.h>\r
+#include "AliExternalTrackParam.h"\r
+#include "AliLog.h"\r
+\r
+//-------------------------------------------------------------------------\r
+// Current support and development: Ruben Shahoyan (Ruben.Shahoyan@cern.ch) \r
+//-------------------------------------------------------------------------\r
+\r
+\r
+class KMCLayer;\r
+class KMCCluster;\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+class KMCProbe : public AliExternalTrackParam {\r
+ public:\r
+  enum {kBitKilled=BIT(14)};\r
+  enum {kNDOF=5};\r
+  enum {kY2=0,kZ2=2,kSnp2=5,kTgl2=9,kPtI2=14};\r
+  enum {kY,kZ,kSnp,kTgl,kPtI};\r
+  //\r
+ KMCProbe() : fMass(0.14),fChi2(0),fHits(0),fFakes(0),fNHits(0),fNHitsITS(0),fNHitsITSFake(0),fInnLrCheck(fgNITSLayers) {for (int i=kMaxITSLr;i--;) fClID[i] =-2;\r
+}\r
+  KMCProbe(KMCProbe& src);\r
+  KMCProbe& operator=(const KMCProbe& src);\r
+  virtual ~KMCProbe() {}\r
+  virtual   void   Print(Option_t* option = "") const;\r
+  virtual   void   Reset() {fMass=0.14; fChi2=0; fHits=fFakes=0;  for (int i=kMaxITSLr;i--;) fClID[i]=-2; AliExternalTrackParam::Reset();}\r
+  virtual   Bool_t IsSortable()                     const {return kTRUE;}\r
+  virtual   Int_t  Compare(const TObject* obj)      const;\r
+  void      ResetCovMat();\r
+  //\r
+  void      Kill(Bool_t v=kTRUE)                          {SetBit(kBitKilled,v);}\r
+  Bool_t    IsKilled()                              const {return TestBit(kBitKilled);}\r
+  //\r
+  Bool_t    CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward=kTRUE);\r
+  Bool_t    GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir=0) const;\r
+  Bool_t    PropagateToR(double r, double b, int dir=0);\r
+  \r
+  //\r
+  void      SetMass(double m=0.14)                      {fMass = m;}\r
+  Double_t  GetMass()                             const {return fMass;}\r
+  Double_t  GetChi2()                             const {return fChi2;}\r
+  void      SetChi2(double chi2)                        {fChi2 = chi2;}\r
+  void      AddChi2(double chi2)                        {fChi2 += chi2;}\r
+  void      SetInnerLrChecked(Int_t n)                  {if (n<fgNITSLayers) fInnLrCheck = n;}\r
+  UShort_t  GetNITSHits()                         const {return fNHitsITS;}\r
+  UShort_t  GetNFakeITSHits()                     const {return fNHitsITSFake;}\r
+  UShort_t  GetNHits()                            const {return fNHits;}\r
+  UShort_t  GetInnerLayerChecked()                const {return fInnLrCheck;}\r
+  UInt_t&   GetHitsPatt()                               {return fHits;}\r
+  UInt_t&   GetFakesPatt()                              {return fFakes;}\r
+  Double_t  GetNormChi2(Bool_t penalize=kFALSE)   const;\r
+  void      AddHit(Int_t lr, double chi2, Int_t clID=-1);\r
+  void      ResetHit(Int_t lr);\r
+  Bool_t    IsHit(Int_t lr)                       const {return (lr<fgNITSLayers) ? IsWBit(fHits,lr)  : kFALSE;}\r
+  Bool_t    IsHitFake(Int_t lr)                   const {return (lr<fgNITSLayers) ? IsWBit(fFakes,lr) : kFALSE;}\r
+  Bool_t    PropagateToCluster(KMCCluster* cl, Double_t b);\r
+  // protected: \r
+  static void   SetWBit(UInt_t &patt,UInt_t bit)               {patt |= 0x1<<bit;}\r
+  static void   ResetWBit(UInt_t &patt,UInt_t bit)             {patt &= ~(0x1<<bit);}\r
+  static Bool_t IsWBit(const UInt_t &patt,const UInt_t bit)    {return patt&(0x1<<bit);}\r
+  static void   SetNITSLayers(Int_t n)                         {fgNITSLayers = n;}\r
+  static int    GetNITSLayers()                                {return fgNITSLayers;}\r
+  //\r
+  static Double_t GetMissingHitPenalty()                        {return fgMissingHitPenalty;}\r
+  static void     SetMissingHitPenalty(double p=2.)             {fgMissingHitPenalty = p;}\r
+  //\r
+ public:\r
+  enum {kMaxITSLr=12};\r
+  Float_t fMass;   // mass\r
+  Float_t fChi2;   // total chi2\r
+  UInt_t  fHits;   // pattern on hits (max 32!)\r
+  UInt_t  fFakes;  // pattern of fakes among hits\r
+  UShort_t fNHits;    // total hits\r
+  UShort_t fNHitsITS; // total ITS hits\r
+  UShort_t fNHitsITSFake; // number of fake ITS hits\r
+  UShort_t fInnLrCheck;   // lowest active layer where update was checked\r
+  Int_t    fClID[kMaxITSLr];\r
+  //\r
+  static Int_t    fgNITSLayers;\r
+  static Double_t fgMissingHitPenalty;\r
+  ClassDef(KMCProbe,2);  \r
+};\r
+\r
+//_______________________________________\r
+inline Double_t KMCProbe::GetNormChi2(Bool_t penalize) const\r
+{\r
+  // normalized chi2, penilized for missing hits\r
+  if (fNHitsITS<3) return 0;\r
+  double chi2 = fChi2;\r
+  if (penalize) {\r
+    int nMiss = fgNITSLayers - fNHitsITS - fInnLrCheck;\r
+    chi2 = fChi2 + fgMissingHitPenalty*nMiss;\r
+  }\r
+  return chi2/( (fNHitsITS<<1)-kNDOF);\r
+}\r
+\r
+//_______________________________________\r
+inline void KMCProbe::AddHit(Int_t lr, double chi2, Int_t clID) {\r
+  // note: lr is active layer ID\r
+  if (lr<0) return;\r
+  fNHits++;\r
+  if (lr<fgNITSLayers) {\r
+    fChi2 += chi2;\r
+    SetWBit(fHits,lr); \r
+    fNHitsITS++;\r
+    if (clID>-1) {\r
+      SetWBit(fFakes,lr);\r
+      fNHitsITSFake++;\r
+    }\r
+    fClID[lr] = clID;\r
+    //else ResetWBit(fFakes,lr);\r
+  }\r
+}\r
+\r
+//_______________________________________\r
+inline void KMCProbe::ResetHit(Int_t lr) {\r
+  // note: lr is active layer ID\r
+  if (lr>=fgNITSLayers) return; \r
+  if (IsWBit(fHits,lr))  {fNHitsITS--;     ResetWBit(fHits,lr);}\r
+  if (IsWBit(fFakes,lr)) {fNHitsITSFake--; ResetWBit(fFakes,lr);}\r
+}\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+class KMCCluster : public TObject {\r
+ public:\r
+  //\r
+  enum {kBitKilled=BIT(14)};\r
+ KMCCluster(Float_t y=0, Float_t z=0, Float_t x=0, Float_t phi=0) : fY(y),fZ(z),fX(x),fPhi(phi) {}\r
+  KMCCluster(KMCCluster &src);\r
+  KMCCluster& operator=(const KMCCluster& src);\r
+  virtual ~KMCCluster() {}\r
+  void Reset()               {Clear();}\r
+  //\r
+  Double_t GetY()    const   {return fY;}\r
+  Double_t GetX()    const   {return fX;}\r
+  Double_t GetZ()    const   {return fZ;}\r
+  Double_t GetPhi()  const   {return fPhi;}\r
+  //\r
+  void    Kill(Bool_t v=kTRUE)          {SetBit(kBitKilled,v);}\r
+  Bool_t  IsKilled()              const {return TestBit(kBitKilled);}\r
+  Float_t fY; \r
+  Float_t fZ; \r
+  Float_t fX;\r
+  Float_t fPhi;\r
+  void Set(Float_t y, Float_t z, Float_t x, Float_t phi) {fY=y; fZ=z; fX=x; fPhi=phi; ResetBit(kBitKilled);}\r
+  virtual void Print(Option_t * = 0) const;\r
+  ClassDef(KMCCluster,1);\r
+};\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+class KMCLayer : public TNamed {\r
+public:\r
+  enum {kTypeNA=-1,kITS,kTPC,kBitVertex=BIT(15)};\r
+  KMCLayer(char *name);\r
+  Float_t GetRadius()   const {return fR;}\r
+  Float_t GetRadL()     const {return fx2X0;}\r
+  Float_t GetXTimesRho() const {return fXRho;}\r
+  Float_t GetPhiRes()   const {return fPhiRes;}\r
+  Float_t GetZRes()     const {return fZRes;}\r
+  Float_t GetLayerEff() const {return fEff;}\r
+  Int_t   GetActiveID() const {return fActiveID;}\r
+  virtual void  Print(Option_t* option = "") const;\r
+  //\r
+  Bool_t IsDead()       const {return fIsDead;}\r
+  Bool_t IsITS()        const {return fType==kITS;}\r
+  Bool_t IsTPC()        const {return fType==kTPC;}\r
+  Bool_t IsVertex()     const {return TestBit(kBitVertex);}\r
+  //\r
+  Int_t    AddBgCluster(double y,double z,double x,double phi) {int n=GetNBgClusters(); new (fClBg[n]) KMCCluster(y,z,x,phi); return ++n;}\r
+  KMCCluster* GetBgCluster(Int_t i) {return (KMCCluster*)fClBg[i];}\r
+  KMCCluster* GetMCCluster()        {return (KMCCluster*)&fClMC;}\r
+  //\r
+  KMCProbe*       GetAnProbe()    const {return (KMCProbe*)&fTrCorr;}\r
+  Int_t           GetNMCTracks()  const {return fTrMC.GetEntries();}\r
+  KMCProbe*       AddMCTrack(KMCProbe* src=0);\r
+  KMCProbe*       GetMCTrack(Int_t it) const {return (KMCProbe*)fTrMC[it];}\r
+  KMCProbe*       GetWinnerMCTrack()  {if (!fTrMC.IsSorted()) fTrMC.Sort(); return fTrMC.GetEntries() ? (KMCProbe*)fTrMC[0]:0;}\r
+  TClonesArray*         GetMCTracks()        const {return (TClonesArray*)&fTrMC;}\r
+  void Reset();\r
+  void ResetMC() { fClBg.Clear(); fTrMC.Clear();}\r
+  void ResetBgClusters() { fClBg.Clear(); }\r
+  void ResetMCTracks()   { fTrMC.Clear(); }\r
+  Int_t GetNBgClusters() const { return fClBg.GetEntries();}\r
+  static Double_t GetDefEff()   {return fgDefEff;}\r
+  static void     SetDefEff(double eff=1) {fgDefEff = eff>1. ? 1.: (eff<0? 0:eff);}\r
+  //\r
+  //\r
+  Float_t fR; \r
+  Float_t fx2X0;\r
+  Float_t fXRho;    // x*density\r
+  Float_t fPhiRes; \r
+  Float_t fZRes;   \r
+  Float_t fEff;\r
+  Bool_t  fIsDead;\r
+  Int_t   fType;   // its, tpc etc\r
+  Int_t   fActiveID;   // active layer id\r
+  Float_t fSig2EstD;\r
+  Float_t fSig2EstZ;\r
+  //\r
+  KMCCluster   fClCorr;     // ideal cluster\r
+  KMCCluster   fClMC;       // MC cluster (from MS scattered track)\r
+  TClonesArray fClBg; // bg clusters for MC\r
+  //\r
+  KMCProbe     fTrCorr;   // ideal track\r
+  TClonesArray fTrMC;      // MC tracks\r
+  //\r
+  static Double_t fgDefEff;\r
+  ClassDef(KMCLayer,1);\r
+};\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+class KMCDetector : public TNamed {\r
+ public:\r
+  enum {kUtilHisto=BIT(14)};\r
+  KMCDetector();\r
+  KMCDetector(char *name,char *title);\r
+  virtual ~KMCDetector();\r
+\r
+  void AddLayer(char *name, Float_t radius, Float_t radL, Float_t xrho=0, Float_t phiRes=999999, Float_t zRes=999999, Float_t eff=-1);\r
+  Int_t GetLayerID(Int_t actID) const;\r
+  void KillLayer(char *name);\r
+  void SetRadius(char *name, Float_t radius);\r
+  void SetRadiationLength(char *name, Float_t radL);\r
+  void SetResolution(char *name, Float_t phiRes=999999, Float_t zRes=999999);\r
+  void SetLayerEfficiency(char *name, Float_t eff=1.0);\r
+  void RemoveLayer(char *name);\r
+\r
+  Float_t GetRadius(char *name);\r
+  Float_t GetRadiationLength(char *name);\r
+  Float_t GetResolution(char *name, Int_t axis=0);\r
+  Float_t GetLayerEfficiency(char *name);\r
+\r
+  void PrintLayout(); \r
+  void PlotLayout(Int_t plotDead = kTRUE);\r
+  \r
+  void MakeAliceAllNew(Bool_t flagTPC =0,  Bool_t flagMon=1,int setVer=0);\r
+  void MakeAliceCurrent(Bool_t flagTPC =0, Int_t AlignResiduals = 0);\r
+  void AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip=1);\r
+  void RemoveTPC();\r
+\r
+  void SetBField(Float_t bfield) {fBFieldG = bfield*10; }\r
+  Float_t GetBField() const {return fBFieldG/10; }\r
+  void SetLhcUPCscale(Float_t lhcUPCscale) {fLhcUPCscale = lhcUPCscale; }\r
+  Float_t GetLhcUPCscale() const { return fLhcUPCscale; }\r
+  void SetParticleMass(Float_t particleMass) {fParticleMass = particleMass; }\r
+  Float_t GetParticleMass() const { return fParticleMass; }\r
+  void SetIntegrationTime(Float_t integrationTime) {fIntegrationTime = integrationTime; }\r
+  Float_t GetIntegrationTime() const { return fIntegrationTime; }\r
+  void SetMaxRadiusOfSlowDetectors(Float_t maxRadiusSlowDet) {fMaxRadiusSlowDet =  maxRadiusSlowDet; }\r
+  Float_t GetMaxRadiusOfSlowDetectors() const { return fMaxRadiusSlowDet; }\r
+  void SetAvgRapidity(Float_t avgRapidity) {fAvgRapidity = avgRapidity; }\r
+  Float_t GetAvgRapidity() const { return fAvgRapidity; }\r
+  void SetConfidenceLevel(Float_t confLevel) {fConfLevel = confLevel; }\r
+  Float_t GetConfidenceLevel() const { return fConfLevel; }\r
+  void SetAtLeastCorr(Int_t atLeastCorr ) {fAtLeastCorr = atLeastCorr; }\r
+  Float_t GetAtLeastCorr() const { return fAtLeastCorr; }\r
+  void SetAtLeastFake(Int_t atLeastFake ) {fAtLeastFake = atLeastFake; }\r
+  Float_t GetAtLeastFake() const { return fAtLeastFake; }\r
+\r
+  void SetdNdEtaCent(Int_t dNdEtaCent ) {fdNdEtaCent = dNdEtaCent; }\r
+  Float_t GetdNdEtaCent() const { return fdNdEtaCent; }\r
+\r
+  Int_t GetNLayers()          const {return fLayers.GetEntries(); }\r
+  Int_t GetNActiveLayers()    const {return fNActiveLayers; }\r
+  Int_t GetNActiveITSLayers() const {return fNActiveITSLayers; }\r
+\r
+  Bool_t SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta);\r
+  Bool_t SolveSingleTrackViaKalmanMC(int offset=6);\r
+  Bool_t SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr=0, int nMC=10000,int offset=6);\r
+  KMCProbe* KalmanSmooth(int actLr, int actMin,int actMax) const;\r
+  KMCProbe* KalmanSmoothFull(int actLr, int actMin,int actMax) const; //TBD\r
+  void   EliminateUnrelated();\r
+  //\r
+  KMCProbe* PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi=0,double x=0,double y=0,double z=0);\r
+  Bool_t TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const;\r
+  void   ApplyMS(KMCProbe* trc,  double x2x0) const;\r
+  Bool_t PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const;\r
+  Bool_t UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster=kTRUE) const;\r
+\r
+  // Helper functions\r
+  Double_t ThetaMCS                 ( Double_t mass, Double_t RadLength, Double_t momentum ) const;\r
+  Double_t ProbGoodHit              ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbGoodChiSqHit         ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
+  Double_t ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )   ; \r
\r
+  // Howard W. hit distribution and convolution integral\r
+  Double_t Dist              ( Double_t Z, Double_t radius ) ;  \r
+  Double_t HitDensity        ( Double_t radius )   ;\r
+  Double_t UpcHitDensity     ( Double_t radius )   ;\r
+  Double_t IntegratedHitDensity  ( Double_t multiplicity, Double_t radius )   ;\r
+  Double_t OneEventHitDensity    ( Double_t multiplicity, Double_t radius ) const   ;\r
+  Double_t D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) const ;\r
+  \r
+  TGraph* GetGraphMomentumResolution(Int_t color, Int_t linewidth=1);\r
+  TGraph* GetGraphPointingResolution(Int_t axis,Int_t color, Int_t linewidth=1);\r
+  TGraph* GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth=1);\r
+\r
+  TGraph* GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth=1);\r
+\r
+  TGraph* GetGraphRecoEfficiency(Int_t particle, Int_t color, Int_t linewidth=1); \r
+  TGraph* GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth);\r
+\r
+  TGraph* GetGraph(Int_t number, Int_t color, Int_t linewidth=1);\r
+\r
+  void MakeStandardPlots(Bool_t add =0, Int_t color=1, Int_t linewidth=1,Bool_t onlyPionEff=0);\r
+\r
+  // method to extend AliExternalTrackParam functionality\r
+  Bool_t IsZero(double val, double tol=1e-9) const {return TMath::Abs(val)<tol;}\r
+  TList *GetLayers()                   const {return (TList*)&fLayers;}\r
+  KMCLayer* GetLayer(Int_t i)          const {return (KMCLayer*) fLayers.At(i);}\r
+  KMCLayer* GetActiveLayer(Int_t actID)    const {int pid=GetLayerID(actID); return pid<0 ? 0:GetLayer(pid);}\r
+  KMCLayer* GetLayer(const char* name) const {return (KMCLayer*) fLayers.FindObject(name);}\r
+  KMCProbe* GetProbeTrack()       const {return (KMCProbe*)&fProbe;}\r
+  void   ClassifyLayers();\r
+  void   ResetMCTracks(int maxLr=-1);                      \r
+  //\r
+  Bool_t   GetUseBackground()               const {return fUseBackground;}\r
+  void     SetUseBackground(Bool_t v=kTRUE)       {fUseBackground = v;}\r
+  void     CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP);\r
+  void     ResetSearchLimits() {fBgYMin=fBgZMin=1e6; fBgYMax=fBgZMax=-1e6; fNBgLimits=0;}\r
+  void     UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr);\r
+  Int_t    GenBgClusters(KMCLayer* lr);\r
+  Bool_t   NeedToKill(KMCProbe* probe) const;\r
+  Double_t PropagateBack(KMCProbe* trc);\r
+  //\r
+  // definition of reconstructable track\r
+  void     RequireMaxChi2Cl(double cut=25.)           {fMaxChi2Cl = cut>0 ? cut:9; fMaxChi2ClSQ = TMath::Sqrt(fMaxChi2Cl);}\r
+  void     RequireMinITSHits(Int_t n=4)               {fMinITSHits = n;}\r
+  void     RequireMaxNormChi2NDF(double cut=5.)       {fMaxNormChi2NDF = cut>0 ? cut:9;}\r
+  void     RequirePattern(UInt_t *patt, int groups);\r
+  //\r
+  Double_t GetMaxChi2Cl()                      const {return fMaxChi2Cl;}\r
+  Double_t GetMaxNormChi2NDFusterKMC()              const {return fMaxNormChi2NDF;}\r
+  Int_t    GetMinITSHits()                     const {return fMinITSHits;}\r
+  //\r
+  Double_t GetUpdCalls()                       const {return fUpdCalls;}\r
+  void     InitMCWatchHistos();\r
+  TH2F*    GetHMCLrResidRPhi()                 const {return fHMCLrResidRPhi;}\r
+  TH2F*    GetHMCLrResidZ()                    const {return fHMCLrResidZ;}\r
+  TH2F*    GetHMCLrChi2()                      const {return fHMCLrChi2;}\r
+  //\r
+  void     PrintITS(Option_t* opt="") const {for (int i=0;i<=fLastActiveITSLayer;i++) if (!GetLayer(i)->IsDead()) GetLayer(i)->Print(opt);}\r
+  static void SetVtxConstraint(double d=-1, double z=-1) {fgVtxConstraint[0]=d; fgVtxConstraint[1]=z;}\r
+  //\r
+  void CalcHardSearchLimits(double dzv);\r
+  void SetMaxSeedToPropagate(Int_t n=300) {fMaxSeedToPropagate = n;}\r
+\r
+ protected:\r
\r
+  Int_t fNLayers;        // total number of layers in the model\r
+  Int_t fNActiveLayers;  // number of active layers in the model\r
+  Int_t fNActiveITSLayers;  // number of active ITS layers in the model\r
+  Int_t fLastActiveLayer;       // id of last active layer\r
+  Int_t fLastActiveITSLayer;    // id of last active ITS layer\r
+  Int_t fLastActiveLayerTracked;    // id of last active layer really used for tracking of given pt\r
+  TList fLayers;                // List of layer pointers\r
+  Float_t fBFieldG;             // Magnetic Field in Gauss (set in Tesla)\r
+  Float_t fLhcUPCscale;         // UltraPeripheralElectrons: scale from RHIC to LHC\r
+  Float_t fIntegrationTime;     // electronics integration time\r
+  Float_t fConfLevel;           // Confidence Level for the tracking\r
+  Float_t fAvgRapidity;         // rapidity of the track (= mean)\r
+  Float_t fParticleMass;        // Particle used for tracking. Standard: mass of pion\r
+  Double_t fMaxRadiusSlowDet;   // Maximum radius for slow detectors.  Fast detectors \r
+                                // and only fast detectors reside outside this radius.\r
+  Int_t fAtLeastCorr;     // min. number of correct hits for the track to be "good"\r
+  Int_t fAtLeastFake;     // min. number of fake hits for the track to be "fake"\r
+\r
+  Int_t fdNdEtaCent;       // Multiplicity\r
+  //\r
+  // reconstruction settings\r
+  Double_t fMaxChi2Cl;   // max cluster-track chi2 \r
+  Double_t fMaxNormChi2NDF;// max chi2/NDF to accept\r
+  Int_t    fMinITSHits;  // min ITS hits in track to accept\r
+  //\r
+  Double_t fMaxChi2ClSQ; // precalculated sqrt(chi2);\r
+  //\r
+  Int_t    fMaxSeedToPropagate; // take 1st fMaxSeedToPropagate seeds from each layer\r
+  // background settings\r
+  Bool_t   fUseBackground; // do we want to simulate background?\r
+  Double_t fBgYMin;      // min Y for current layer bg generation\r
+  Double_t fBgYMax;      // max Y for current layer bg generation\r
+  Double_t fBgZMin;      // min Z for current layer bg generation\r
+  Double_t fBgZMax;      // max Z for current layer bg generation\r
+  TArrayD  fBgYMinTr;    // min Y for each seed at current layer\r
+  TArrayD  fBgYMaxTr;    // max Y for each seed at current layer\r
+  TArrayD  fBgZMinTr;    // min Z for each seed at current layer\r
+  TArrayD  fBgZMaxTr;    // max Z for each seed at current layer\r
+  Int_t    fNBgLimits;   // depth of limits array\r
+  //\r
+  enum {kMaxNDetectors = 200};\r
\r
+  Double_t fTransMomenta[40];                          // array of transverse momenta\r
+  Double_t fMomentumRes[40];                           // array of momentum resolution\r
+  Double_t fResolutionRPhi[40];                        // array of rphi resolution\r
+  Double_t fResolutionZ[40];                           // array of z resolution\r
+  Double_t fDetPointRes[kMaxNDetectors][40];           // array of rphi resolution per layer\r
+  Double_t fDetPointZRes[kMaxNDetectors][40];          // array of z resolution per layer\r
+  Double_t fEfficiency[3][40];                         // efficiency for different particles\r
+  Double_t fFake[3][40];                               // fake prob for different particles\r
+  //\r
+  KMCProbe fProbe;\r
+  Double_t fDensFactorEta;                             // density scaling for non-0 eta (precalculated\r
+  //\r
+  Double_t fUpdCalls;                  // number of kalman updates\r
+  TH2F*  fHMCLrResidRPhi;              // Residials on lr, rphi\r
+  TH2F*  fHMCLrResidZ;                 // Residials on lr, rphi\r
+  TH2F*  fHMCLrChi2;                   // track to clusyer chi2 on each lr, rphi\r
+  TArrayI fPattITS;                     // bit pattern of each group of active layers where hit is required\r
+  //\r
+  static Double_t fgVtxConstraint[2];  // if both positive, the vertex is used as constraint (accounted in chi2 but not in update)\r
+  ClassDef(KMCDetector,1);\r
+};\r
+\r
+//____________________________________________________________________________\r
+inline Bool_t KMCProbe::PropagateToCluster(KMCCluster* cl, double b)\r
+{\r
+  // propagate track to cluster frame\r
+  if (!Rotate(cl->GetPhi()) || !PropagateTo(cl->GetX(),b)) {\r
+    AliDebug(2,Form("Failed to propager track to cluster at phi=%.3f X=%.3f",cl->GetPhi(),cl->GetX()));\r
+    if (AliLog::GetGlobalDebugLevel()>1) Print();\r
+    return kFALSE;\r
+  }\r
+  return kTRUE;\r
+}\r
+\r
+//////////////////////////////////////////////////////////////////////////////////////////////\r
+//--------------------------------------------------------------------------------------------\r
+//\r
+// Class to collect summary of certrain track performance\r
+// The track can be selected according to min and max ITS clusters to accept (total, fake and or correct)\r
+// and according to presence of (any and/or correct)  or absence of fakes clusters in certain groups\r
+// of ITS layers.\r
+\r
+class KMCTrackSummary : public TNamed\r
+{\r
+ public:\r
+  KMCTrackSummary(const char* name=0,const char* title=0,int nlr=0);\r
+  virtual ~KMCTrackSummary();\r
+  //\r
+  void Add(const KMCTrackSummary* src, double scl=1);\r
+\r
+  virtual void  Print(Option_t* option = "") const;\r
+  KMCTrackSummary* MakeCopy(const char* pref) const;\r
+  void   SetNamePrefix(const char* pref="");\r
+\r
+  Bool_t CheckTrack(KMCProbe* trc);\r
+  void   AddTrack(KMCProbe* trc);\r
+  //\r
+  TH1*   GetHMCChi2()               const {return fHMCChi2;}\r
+  TH1*   GetHMCSigDCARPhi()         const {return fHMCSigDCARPhi;}\r
+  TH1*   GetHMCSigDCAZ()            const {return fHMCSigDCAZ;}\r
+  TH1*   GetHMCSigPt()              const {return fHMCSigPt;}\r
+  TH2*   GetHMCNCl()                const {return fHMCNCl;}\r
+  TH1*   GetHMCFakePatt()           const {return fHMCFakePatt;}\r
+  //\r
+  Double_t GetCountAcc()            const {return fCountAcc;}\r
+  Double_t GetCountTot()            const {return fCountTot;}\r
+  Double_t GetEff()                 const {return fCountTot>0 ? fCountAcc/fCountTot : 0.;}\r
+  Double_t GetEffErr()              const;\r
+  KMCProbe& GetAnProbe()      const {return (KMCProbe&)fAnProbe;}\r
+  KMCProbe& GetRefProbe()     const {return (KMCProbe&)fRefProbe;}\r
+  void SetAnProbe(KMCProbe* pr)     {if (pr) fAnProbe = *pr;}\r
+  void SetRefProbe(KMCProbe* pr)    {if (pr) fRefProbe = *pr;}\r
+  //\r
+  void SetMinMaxClITS(Int_t mn=0,Int_t mx=999)     {fMinMaxClITS[0]=mn; fMinMaxClITS[1]=mx;}\r
+  void SetMinMaxClITSFake(Int_t mn=0,Int_t mx=999) {fMinMaxClITSFake[0]=mn; fMinMaxClITSFake[1]=mx;}\r
+  void SetMinMaxClITSCorr(Int_t mn=0,Int_t mx=999) {fMinMaxClITSCorr[0]=mn; fMinMaxClITSCorr[1]=mx;}\r
+  //\r
+  Double_t GetUpdCalls()           const {return fCountTot>0 ? fUpdCalls/fCountTot : 0;}\r
+  void     AddUpdCalls(double v)         {fUpdCalls += v;}\r
+  //\r
+  void AddPatternITS(Int_t patt)           {AddPattern(patt,fPattITS);}\r
+  void AddPatternITSFakeExcl(Int_t patt)   {AddPattern(patt,fPattITSFakeExcl);}\r
+  void AddPatternITSCorr(Int_t patt)       {AddPattern(patt,fPattITSCorr);}\r
+  static UInt_t Bits(Bool_t l0=0, Bool_t l1=0, Bool_t l2=0, Bool_t l3=0, Bool_t l4=0, Bool_t l5=0, Bool_t l6=0, Bool_t l7=0,Bool_t  l8=0, Bool_t l9=0,\r
+                    Bool_t l10=0,Bool_t l11=0,Bool_t l12=0,Bool_t l13=0,Bool_t l14=0,Bool_t l15=0,Bool_t l16=0,Bool_t l17=0,Bool_t l18=0,Bool_t l19=0,\r
+                    Bool_t l20=0,Bool_t l21=0,Bool_t l22=0,Bool_t l23=0,Bool_t l24=0,Bool_t l25=0,Bool_t l26=0,Bool_t l27=0,Bool_t l28=0,Bool_t l29=0,\r
+                    Bool_t l30=0,Bool_t l31=0);\r
+  //\r
+  // protected:\r
+  void   AddPattern(Int_t patt,TArrayI& arr) {int n=arr.GetSize(); arr.Set(n+1); arr[n]=patt;}\r
+  Bool_t CheckPattern(UInt_t patt, Int_t cond) {return ((UInt_t)cond)&patt;}\r
+  Int_t  OutOfRange(Int_t n, Int_t range[2])  {if (n<range[0]) return -1; if (n>range[1]) return 1; return 0;}\r
+  static void   PrependTNamed(TNamed* obj, const char* nm=0, const char* tit=0);\r
+  //\r
+ private:\r
+ KMCTrackSummary(const KMCTrackSummary& src) :TNamed(src) {}\r
+ protected:\r
+  //\r
+  Int_t   fNITSLayers;                  // number of its layers\r
+  // track selection conditions\r
+  Int_t   fMinMaxClITS[2];              // min and max N ITS clusters (total)\r
+  Int_t   fMinMaxClITSFake[2];          // min and max N ITS fake clusters  \r
+  Int_t   fMinMaxClITSCorr[2];          // min and max N ITS correct clusters\r
+  //\r
+  TArrayI fPattITS;                     // bit pattern of each group of layers where hit is required\r
+  TArrayI fPattITSFakeExcl;             // bit pattern of each group of layers where fake hit is not tolerated\r
+  TArrayI fPattITSCorr;                 // bit pattern of each group of layers where correc hit is required\r
+  //\r
+  TH1F*  fHMCChi2;                     // track chi2\r
+  TH1F*  fHMCSigDCARPhi;               // DCA distribution in RPhi from MC\r
+  TH1F*  fHMCSigDCAZ;                  // DCA distribution in Z for from MC\r
+  TH1F*  fHMCSigPt;                    // Pt difference distribution from MC\r
+  TH2F*  fHMCNCl;                      // number of clusters\r
+  TH1F*  fHMCFakePatt;                 // Fakes pattern\r
+  //\r
+  Double_t fCountAcc;                  // track counter of accepted tracks\r
+  Double_t fCountTot;                  // track counter of all tracks\r
+  Double_t fUpdCalls;                  // total update calls\r
+  //\r
+  KMCProbe fRefProbe;                  // reference probe\r
+  KMCProbe fAnProbe;                   // analital probe\r
+  static Int_t fgSumCounter;           // global counter of bins\r
+\r
+  ClassDef(KMCTrackSummary,1)\r
+};\r
+\r
+inline Double_t KMCTrackSummary::GetEffErr() const {\r
+  //\r
+  if (fCountTot<1) return 0;\r
+  double r = fCountAcc/fCountTot;\r
+  double err = r*(1.-r)/fCountTot;\r
+  return err>0 ? TMath::Sqrt(err) : 0;\r
+}\r
+\r
+inline KMCTrackSummary* KMCTrackSummary::MakeCopy(const char* pref) const {\r
+  // create a copy, prepending all names with prefix\r
+  KMCTrackSummary* sm = (KMCTrackSummary*) Clone(this->GetName()); \r
+  sm->SetNamePrefix(pref); \r
+  return sm;\r
+} \r
+\r
+#endif\r
diff --git a/ITS/UPGRADE/PrepSummaryKMC.C b/ITS/UPGRADE/PrepSummaryKMC.C
new file mode 100755 (executable)
index 0000000..517f5c0
--- /dev/null
@@ -0,0 +1,195 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TROOT.h>
+#include <TFile.h>
+#include <TDirectory.h>
+#include <TF1.h>
+#include <TH1.h>
+#include <TObjArray.h>
+#include <TGraphErrors.h>
+#include "KMCDetector.h"
+#endif
+
+// Process the output file of testDetKMC.C macro to prepare summary 
+// for the "summary class" icl (definition of acceptable track).
+// See ProcessSummary comment for details.
+// The names of the graphs are prefixed by `pref`.
+// If outF is provided, the graphs will be stored there
+
+
+TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref=0, const char* outF=0);
+TObjArray* ProcessSummary(TObjArray* sums, int icl, const char* pref=0);
+
+enum {kSigAD,kSigAZ,kSigAPt,kSigD,kSigZ,kSigPt,kEff,kUpd};
+TF1* gs = 0;
+
+TObjArray* PrepSummaryKMC(const char* sumf, int cls, const char* pref, const char* outF)
+{
+  if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
+  TFile* fl = TFile::Open(sumf);
+  if (!fl) {printf("No file %s\n",sumf); return 0;}
+  TObjArray* arrs = (TObjArray*)gDirectory->Get("trSum");
+  if (!arrs) {printf("No summary in file %s\n",sumf); return 0;}
+  //
+  TObjArray* sums =  ProcessSummary(arrs,cls,pref);
+  if (!sums) return 0;
+  //
+  if (outF) {
+    TFile* flOut = TFile::Open(outF,"update");
+    if (!flOut) {printf("Failed to open output file %s\n",outF);}
+    else {
+      flOut->WriteObject(sums,sums->GetName(),"kSingleKey");
+      flOut->Close();
+      delete flOut;
+      printf("Stored array %s in %s\n",sums->GetName(), outF);
+    }
+  }
+  //
+  return sums;
+}
+
+
+TObjArray* ProcessSummary(TObjArray* arrs, int icl, const char* pref)
+{
+  // Process TObjArray (e.g. for set of pt bins) of TObjArray of KMCTrackSummary objects:
+  // pick the KMCTrackSummary for "summary class" icl (definition of acceptable track) and create
+  // graphs vs bin.
+  // These graphs are returned in new TObjArray
+  //
+  TString prefs = pref;
+  if (!gs) gs = new TF1("gs","gaus",-1,1);
+  //
+  int nb = arrs->GetEntriesFast();
+  TObjArray* sums = (TObjArray*) arrs->At(0);
+  int nclass = sums->GetEntriesFast();
+  if (icl>=nclass) {printf("summary set has %d classes only, %d requested\n",nclass,icl);return 0;}
+  //
+  KMCTrackSummary* sm = (KMCTrackSummary*)sums->At(icl);
+  //
+  TH1* h;
+  //
+  h = sm->GetHMCSigDCARPhi();  // MC resolution in transverse DCA 
+  TGraphErrors * grSigD = 0;
+  if (h) {
+    grSigD = new TGraphErrors(nb);
+    grSigD->SetName(Form("%s%s",prefs.Data(),h->GetName()));
+    grSigD->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
+  }
+  //
+  TGraphErrors * grSigZ = 0;
+  h = sm->GetHMCSigDCAZ();    // MC resolution in Z DCA
+  if (h) {
+    grSigZ = new TGraphErrors(nb);
+    grSigZ->SetName(Form("%s%s",prefs.Data(),h->GetName()));
+    grSigZ->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
+  }
+  //
+  TGraphErrors * grSigAD = 0; // anaitical estimate for resolution in transverse DCA 
+  {
+    grSigAD = new TGraphErrors(nb);
+    grSigAD->SetName(Form("%s%s",prefs.Data(),"sigmaDan"));
+    grSigAD->SetTitle(Form("%s%s",prefs.Data(),"#sigmaD an"));
+  }
+  //
+  TGraphErrors * grSigAZ = 0; // anaitical estimate for resolution in Z DCA 
+  {
+    grSigAZ = new TGraphErrors(nb);
+    grSigAZ->SetName(Form("%s%s",prefs.Data(),"sigmaZan"));
+    grSigAZ->SetTitle(Form("%s%s",prefs.Data(),"#sigmaZ an"));
+  }
+  //
+  //
+  TGraphErrors * grSigPt = 0; // MC res. in pt
+  {
+    grSigPt = new TGraphErrors(nb);
+    grSigPt->SetName(Form("%s%s",prefs.Data(),"sigmaPt"));
+    grSigPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt"));
+  }
+  //
+  TGraphErrors * grSigAPt = 0; // analitycal res. in pt
+  {
+    grSigAPt = new TGraphErrors(nb);
+    grSigAPt->SetName(Form("%s%s",prefs.Data(),"sigmaPtan"));
+    grSigAPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt an"));
+  }
+
+  //
+  TGraphErrors * grEff = 0; // MC efficiency
+  {
+    grEff = new TGraphErrors(nb);
+    grEff->SetName(Form("%s_rate",prefs.Data()));
+    grEff->SetTitle(Form("%s Rate",prefs.Data()));
+  }
+  //
+  TGraphErrors * grUpd = 0; // number of Kalman track updates
+  {
+    grUpd = new TGraphErrors(nb);
+    grUpd->SetName(Form("%s_updCalls",prefs.Data()));
+    grUpd->SetTitle(Form("%s Updates",prefs.Data()));
+  }
+  //
+  for (int ib=0;ib<nb;ib++) {
+    sums = (TObjArray*) arrs->At(ib);
+    sm = (KMCTrackSummary*)sums->At(icl);
+    KMCProbe& prbRef = sm->GetRefProbe();
+    KMCProbe& prbAn  = sm->GetAnProbe();
+  
+    double pt = prbRef.Pt();
+    //
+    if (grSigAD) {
+      grSigAD->SetPoint(ib, pt,prbAn.GetSigmaY2()>0 ? TMath::Sqrt(prbAn.GetSigmaY2()) : 0.);
+    }
+    //
+    if (grSigAZ) {
+      grSigAZ->SetPoint(ib, pt,prbAn.GetSigmaZ2()>0 ? TMath::Sqrt(prbAn.GetSigmaZ2()) : 0.);
+    }
+    //
+    if (grSigAPt) {
+      double pts = TMath::Sqrt(prbAn.GetSigma1Pt2());
+      grSigAPt->SetPoint(ib, pt,pts>0 ? pts*pt : 0.);
+    }
+    //
+    if (grSigPt) {
+      h = sm->GetHMCSigPt();
+      h->Fit(gs,"0q");
+      grSigPt->SetPoint(ib, pt, gs->GetParameter(2));
+      grSigPt->SetPointError(ib, 0, gs->GetParError(2));
+    }
+    //
+     if (grSigD) {
+      h = sm->GetHMCSigDCARPhi();
+      h->Fit(gs,"0q");
+      grSigD->SetPoint(ib, pt,gs->GetParameter(2));
+      grSigD->SetPointError(ib, 0,gs->GetParError(2));      
+    }
+    //
+    if (grSigZ) {
+      h = sm->GetHMCSigDCAZ();
+      h->Fit(gs,"0q");
+      grSigZ->SetPoint(ib, pt,gs->GetParameter(2));
+      grSigZ->SetPointError(ib, 0,gs->GetParError(2));      
+    }
+    //
+    if (grEff) {
+      grEff->SetPoint(ib, pt,sm->GetEff());
+      grEff->SetPointError(ib, 0,sm->GetEffErr());
+    }
+    //
+    if (grUpd) {
+      grUpd->SetPoint(ib, pt,sm->GetUpdCalls());
+      grUpd->SetPointError(ib, 0, 0);
+    }
+  }
+  //
+  TObjArray* dest = new TObjArray();
+  dest->AddAtAndExpand(grSigAD,kSigAD);
+  dest->AddAtAndExpand(grSigAZ,kSigAZ);  
+  dest->AddAtAndExpand(grSigAPt,kSigAPt);  
+  dest->AddAtAndExpand(grSigD,kSigD);
+  dest->AddAtAndExpand(grSigZ,kSigZ);  
+  dest->AddAtAndExpand(grSigPt,kSigPt);  
+  dest->AddAtAndExpand(grEff,kEff);
+  dest->AddAtAndExpand(grUpd,kUpd);
+  //
+  if (!prefs.IsNull()) dest->SetName(pref);
+  return dest;
+}
diff --git a/ITS/UPGRADE/prepPlots1KMC.C b/ITS/UPGRADE/prepPlots1KMC.C
new file mode 100755 (executable)
index 0000000..1d32303
--- /dev/null
@@ -0,0 +1,303 @@
+{
+  gROOT->LoadMacro("KMCDetector.cxx+");
+  gROOT->LoadMacro("PrepSummaryKMC.C+");
+  //
+  TObjArray* NewSt100Perf       = PrepSummaryKMC("NewStEff100NoVtx.root",0,"NewStEff100Rerf");
+  TObjArray* NewSt100Cor5       = PrepSummaryKMC("NewStEff100NoVtx.root",1,"NewStEff100Cor5");
+  TObjArray* NewSt100Cor5hs2in  = PrepSummaryKMC("NewStEff100NoVtx.root",2,"NewStEff100Cor5hs2in");
+  TObjArray* NewSt100Cor5hd2in  = PrepSummaryKMC("NewStEff100NoVtx.root",3,"NewStEff100Cor5hd2in");
+  TObjArray* NewSt100Cor5hs2inG = PrepSummaryKMC("NewStEff100NoVtx.root",4,"NewStEff100Cor5hs2inG");
+  TObjArray* NewSt100Cor5hd2inG = PrepSummaryKMC("NewStEff100NoVtx.root",5,"NewStEff100Cor5hd2inG");
+  TObjArray* NewSt100Fake       = PrepSummaryKMC("NewStEff100NoVtx.root",6,"NewStEff100Fake");
+
+  TObjArray* NewSt095Perf       = PrepSummaryKMC("NewStEff095NoVtx.root",0,"NewStEff095Rerf");
+  TObjArray* NewSt095Cor5       = PrepSummaryKMC("NewStEff095NoVtx.root",1,"NewStEff095Cor5");
+  TObjArray* NewSt095Cor5hs2in  = PrepSummaryKMC("NewStEff095NoVtx.root",2,"NewStEff095Cor5hs2in");
+  TObjArray* NewSt095Cor5hd2in  = PrepSummaryKMC("NewStEff095NoVtx.root",3,"NewStEff095Cor5hd2in");
+  TObjArray* NewSt095Cor5hs2inG = PrepSummaryKMC("NewStEff095NoVtx.root",4,"NewStEff095Cor5hs2inG");
+  TObjArray* NewSt095Cor5hd2inG = PrepSummaryKMC("NewStEff095NoVtx.root",5,"NewStEff095Cor5hd2inG");
+  TObjArray* NewSt095Fake       = PrepSummaryKMC("NewStEff095NoVtx.root",6,"NewStEff095Fake");
+
+  //
+  TObjArray* NewRS100Perf       = PrepSummaryKMC("NewRSEff100NoVtx.root",0,"NewRSEff100Rerf");
+  TObjArray* NewRS100Cor5       = PrepSummaryKMC("NewRSEff100NoVtx.root",1,"NewRSEff100Cor5");
+  TObjArray* NewRS100Cor5hs2in  = PrepSummaryKMC("NewRSEff100NoVtx.root",2,"NewRSEff100Cor5hs2in");
+  TObjArray* NewRS100Cor5hd2in  = PrepSummaryKMC("NewRSEff100NoVtx.root",3,"NewRSEff100Cor5hd2in");
+  TObjArray* NewRS100Cor5hs2inG = PrepSummaryKMC("NewRSEff100NoVtx.root",4,"NewRSEff100Cor5hs2inG");
+  TObjArray* NewRS100Cor5hd2inG = PrepSummaryKMC("NewRSEff100NoVtx.root",5,"NewRSEff100Cor5hd2inG");
+  TObjArray* NewRS100Fake       = PrepSummaryKMC("NewRSEff100NoVtx.root",6,"NewRSEff100Fake");
+
+  TObjArray* NewRS095Perf       = PrepSummaryKMC("NewRSEff095NoVtx.root",0,"NewRSEff095Rerf");
+  TObjArray* NewRS095Cor5       = PrepSummaryKMC("NewRSEff095NoVtx.root",1,"NewRSEff095Cor5");
+  TObjArray* NewRS095Cor5hs2in  = PrepSummaryKMC("NewRSEff095NoVtx.root",2,"NewRSEff095Cor5hs2in");
+  TObjArray* NewRS095Cor5hd2in  = PrepSummaryKMC("NewRSEff095NoVtx.root",3,"NewRSEff095Cor5hd2in");
+  TObjArray* NewRS095Cor5hs2inG = PrepSummaryKMC("NewRSEff095NoVtx.root",4,"NewRSEff095Cor5hs2inG");
+  TObjArray* NewRS095Cor5hd2inG = PrepSummaryKMC("NewRSEff095NoVtx.root",5,"NewRSEff095Cor5hd2inG");
+  TObjArray* NewRS095Fake       = PrepSummaryKMC("NewRSEff095NoVtx.root",6,"NewRSEff095Fake");
+
+  gStyle->SetOptStat(0);
+  gStyle->SetOptTitle(0);
+  TCanvas* c1 = new TCanvas();
+  c1->Draw();
+  TH1F* heff = new TH1F("heff","",100,0.1,1.7);
+  heff->GetYaxis()->SetNdivisions(521);
+  heff->SetXTitle("p_{T}");
+  heff->SetYTitle("Rate");  
+  heff->SetMinimum(0);
+  heff->SetMaximum(1.09);  
+  heff->Draw();
+  gPad->SetGrid(1,1);
+  //
+  TGraphErrors* gr;
+  //
+  Int_t col = 0;
+  Int_t mtype = 20;
+  Int_t msize = 1;
+  TLegend *leg1 = new TLegend(0.20,0.25,0.95,0.65," ");
+  leg1->SetNColumns(2);
+  leg1->SetFillColor(0);
+  TLegendEntry* le=0;
+  //
+
+  col = kRed;
+  mtype = 20;
+  gr = (TGraphErrors*)NewSt100Cor5->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 100%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kRed+2;
+  mtype = 21;
+  gr = (TGraphErrors*)NewSt100Cor5hs2in->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 100%, min 5h(2 inner)","lp");
+  le->SetTextColor(col);
+  // 
+  //
+  //-------------------------
+
+  col = kRed;
+  mtype = 24;
+  gr = (TGraphErrors*)NewSt095Cor5->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kRed+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewSt095Cor5hs2in->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 95%, min 5h(2 inner)","lp");
+  le->SetTextColor(col);
+  // 
+  //
+  // 
+  col = kRed;
+  mtype = 20;
+  gr = (TGraphErrors*)NewSt100Fake->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->SetLineStyle(2);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 100%, Fakes","lp");
+  le->SetTextColor(col);
+  //
+  // 
+  col = kRed+2;
+  mtype = 21;
+  gr = (TGraphErrors*)NewSt095Fake->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->SetLineStyle(2);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New St, eff 95%, Fakes","lp");
+  le->SetTextColor(col);
+  //
+  //
+
+  //-------------------------------------------------------------------------------------
+  //
+  
+  col = kBlue;
+  mtype = 20;
+  gr = (TGraphErrors*)NewRS100Cor5->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 100%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kGreen+2;
+  mtype = 21;
+  gr = (TGraphErrors*)NewRS100Cor5hd2in->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 100%, min 5h(2 inner pair)","lp");
+  le->SetTextColor(col);
+  // 
+  //
+  //-------------------------
+
+  col = kBlue;
+  mtype = 24;
+  gr = (TGraphErrors*)NewRS095Cor5->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kGreen+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewRS095Cor5hd2in->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 95%, min 5h(2 inner pair)","lp");
+  le->SetTextColor(col);
+  // 
+  //
+  // 
+  col = kBlue;
+  mtype = 20;
+  gr = (TGraphErrors*)NewRS100Fake->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->SetLineStyle(2);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 100%, Fakes","lp");
+  le->SetTextColor(col);
+  //
+  // 
+  col = kGreen+2;
+  mtype = 21;
+  gr = (TGraphErrors*)NewRS095Fake->At(kEff);
+  SetGStyle(gr,col,mtype,msize);
+  gr->SetLineStyle(2);
+  gr->Draw("pc");
+  le = leg1->AddEntry(gr,"New+2in dbl, eff 95%, Fakes","lp");
+  le->SetTextColor(col);
+  //
+  leg1->Draw();
+  //
+  SaveCanvas(c1,"fig/effFakes_new_st_2indbl","cg");
+
+  //===============================================
+  TCanvas* c2 = new TCanvas();
+  c2->Draw();
+  gPad->Modified();
+  gPad->SetLeftMargin(0.15);
+  gPad->Modified();
+  TLegend *leg2 = new TLegend(0.45,0.5,0.9,0.85,"");
+  leg2->SetFillColor(0);
+  
+  TH1F* hsigD = new TH1F("sigd","",100,0.1,1.7);
+  hsigD->SetXTitle("p_{T}");
+  hsigD->SetYTitle("resolution, #sigma_{r#phi}");  
+  hsigD->SetMinimum(0);
+  hsigD->SetMaximum(0.025);  
+  hsigD->Draw();
+  hsigD->GetYaxis()->SetTitleOffset(1.4);
+  hsigD->GetYaxis()->SetTitleSize(0.05);
+  // 
+  gPad->SetGrid(1,1);
+  //
+  //------------------------------------------------------------
+
+  col = kRed;
+  mtype = 20;
+  gr = (TGraphErrors*)NewSt095Cor5->At(kSigD);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg2->AddEntry(gr,"New St, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kRed+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewSt095Cor5hs2in->At(kSigD);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg2->AddEntry(gr,"New St, eff 95%, min 5h(2 inner)","lp");
+  le->SetTextColor(col);
+  // 
+  //--------------------
+  col = kBlue;
+  mtype = 20;
+  gr = (TGraphErrors*)NewRS095Cor5->At(kSigD);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg2->AddEntry(gr,"New+2in dbl, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kGreen+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewRS095Cor5hd2in->At(kSigD);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg2->AddEntry(gr,"New+2in dbl, eff 95%, min 5h(2 inner pair)","lp");
+  le->SetTextColor(col);
+  // 
+  leg2->Draw();
+  SaveCanvas(c2,"fig/resD_new_st_2indbl","cg");
+
+
+  //----------------------------------------------
+  TCanvas* c3 = new TCanvas();
+  c3->Draw();
+  gPad->Modified();
+  gPad->SetLeftMargin(0.15);
+  gPad->Modified();
+  TLegend *leg3 = new TLegend(0.45,0.5,0.9,0.85,"");
+  leg3->SetFillColor(0);
+  
+  TH1F* hsigZ = new TH1F("sigd","",100,0.1,1.7);
+  hsigZ->SetXTitle("p_{T}");
+  hsigZ->SetYTitle("resolution, #sigma_{Z}");  
+  hsigZ->SetMinimum(0);
+  hsigZ->SetMaximum(0.025);  
+  hsigZ->Draw();
+  hsigZ->GetYaxis()->SetTitleOffset(1.4);
+  hsigZ->GetYaxis()->SetTitleSize(0.05);
+  // 
+  gPad->SetGrid(1,1);
+  //
+
+
+  col = kRed;
+  mtype = 20;
+  gr = (TGraphErrors*)NewSt095Cor5->At(kSigZ);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg3->AddEntry(gr,"New St, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kRed+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewSt095Cor5hs2in->At(kSigZ);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg3->AddEntry(gr,"New St, eff 95%, min 5h(2 inner)","lp");
+  le->SetTextColor(col);
+  // 
+  //--------------------
+  col = kBlue;
+  mtype = 20;
+  gr = (TGraphErrors*)NewRS095Cor5->At(kSigZ);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg3->AddEntry(gr,"New+2in dbl, eff 95%, min 5h","lp");
+  le->SetTextColor(col);
+  //
+  col = kGreen+2;
+  mtype = 25;
+  gr = (TGraphErrors*)NewRS095Cor5hd2in->At(kSigZ);
+  SetGStyle(gr,col,mtype,msize);
+  gr->Draw("pc");
+  le = leg3->AddEntry(gr,"New+2in dbl, eff 95%, min 5h(2 inner pair)","lp");
+  le->SetTextColor(col);
+  // 
+  //
+  leg3->Draw();
+  SaveCanvas(c3,"fig/resZ_new_st_2indbl","cg");
+
+}
diff --git a/ITS/UPGRADE/readmeKMC b/ITS/UPGRADE/readmeKMC
new file mode 100644 (file)
index 0000000..e4f73d7
--- /dev/null
@@ -0,0 +1,44 @@
+Here are some instructions for running:
+Once the KMCDetector object is initialized (similar to FER, see "Detector, tracking initialization"
+part of the testDetKMC.C macro) one should define for each kinematic bin to test a TObjArray with
+KMCTrackSummary objects (done in the testDetKMC.C by the CreateTrackConditions routine):
+this object is basically a set of standard histos (with DCA, pt_rec - pt_gen distibutions etc + stat. of generated
+and found tracks) which are filled provided the reco track satisfies to user defined conditions, like the number of
+hits (one can specify correct or fake hits, specific layers pattern etc). For instance:
+
+sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr); // create a summary object with this name
+
+sm->SetMinMaxClITS(4);             // require at least 4 cluster in the track (the upper limit here is default 999 - no limit)
+
+sm->SetMinMaxClITSFake(0,0); // set min number of fake clusters required (used for fakes studies) and their max number tolerated
+
+sm->SetNamePrefix(nm);          // this histo names will be prefixed by nm
+
+sm->AddPatternITSCorr( sm->Bits(1,1));    // require that the innermost 2 double layers must contribute at least 1 correct cluster
+
+sm->AddPatternITSCorr( sm->Bits(0,0,1));  // require the 3d layer cluster to be correct
+
+sm->AddPatternITSCorr( sm->Bits(0,0,0,1,1)); // measurement in the middle: ask at leas 1 correct cluster in 2 middle layers
+
+sm->AddPatternITSCorr( sm->Bits(0,0,0, 0,0, 1,1)); // same for outer 2 layers
+
+One can define and put in the TObjArray many such objects with different conditions, they will be filled independently.
+The simulation for single kinem. bin is done via call:
+
+its.SolveSingleTrack(mass, pt, eta, trSum, nevR);
+This will generate nevR tracks of given mass, pt, eta, propagate them, reconstruct under imposed bg, sensor err. conditions
+and fill the KMCTrackSummary objects provided via trSum TObjArray (i.e. histos of each object are filled if given reco track
+satisfies the criteria defined for this object).
+
+In the testDetKMC trSum arrays for each pt bin are added to another TObjArray (e.g. summary) which is then stored in
+the sumOut file.
+
+I also attach the PrepSummaryKMC.C macro to fill/return from this output file the TObjArray of TGraphErrors corresponding
+to given id of KMCTrackSummary. The call
+PrepSummaryKMC(const char* sumf, int cls, const char* pref=0, const char* outF=0);
+will open sumf output of the testDetKMC, and  for all pt bins (indexing the entries of the final "summary" TObjArray
+of testDetKMC) it will extract the array of KMCTrackSummary objects defined in the CreateTrackConditions
+for this bin, then extract from it the cls-th KMCTrackSummary object and put the performance values for this pt bin
+(resolutions, eff ) as an entry in the graph.
+
+prepPlotsKMC1.C is just an example of preparing the graphs from the testDetKMC.C output objects.
diff --git a/ITS/UPGRADE/testDetKMC.C b/ITS/UPGRADE/testDetKMC.C
new file mode 100644 (file)
index 0000000..0396ef2
--- /dev/null
@@ -0,0 +1,192 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TCanvas.h>
+#include <TString.h>
+#include <TROOT.h>
+#include <TSystem.h>
+#include <TFile.h>
+#include <TObjArray.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include "KMCDetector.h"
+#endif
+
+TObjArray* CreateTrackConditions(const char* nm);
+TObjArray  summary;
+
+
+TObject* gDet=0;
+
+// ptBins to generate
+double pts[] = {0.15,0.175,0.2,0.225,0.25,0.3,0.4,0.5,0.8,1.1,1.6,2,4,8,13,20};
+//double pts[] = {0.2, 1., 10.};
+
+void testDetKMC(int nev=10000,  // n events to generate per bin
+               int setVer=0,   // optional version to pass for detector building
+               double mass=0.14,  // particle mass to generate
+               double eta=0.45,   // eta to test
+               Bool_t aliceNew=kTRUE,Bool_t tpc=kFALSE, Double_t eff=0.95,  // detector settings
+               Double_t vtxConstr=-1., // use vertex constraint in the fit
+               const char* sumOut="trsumDef.root" // output file name
+               )
+{
+  if (!gROOT->GetClass("KMCDetector")) gROOT->LoadMacro("KMCDetector.cxx+");
+  //
+  // Detector, tracking initialization >>>
+  //
+  if (vtxConstr>0) KMCDetector::SetVtxConstraint(vtxConstr,vtxConstr);
+  KMCLayer::SetDefEff(eff);
+  //
+  KMCDetector *itsdet = new KMCDetector((char*)"ALICE",(char*)"ITS");
+  gDet = itsdet;
+  KMCDetector &its = *itsdet;
+  its.SetAvgRapidity(eta);
+  if (aliceNew)     its.MakeAliceAllNew(tpc,1, setVer); // its sa
+  else              its.MakeAliceCurrent(tpc,0); // its sa  
+  //its.MakeAliceCurrent(0,1); // with tpc
+  its.SetMaxSeedToPropagate(300);
+  its.SetUseBackground();
+  //
+  // min hits per track to validate it
+  its.RequireMinITSHits( its.GetNActiveITSLayers() - 2 );//(int)TMath::Max(4.0,(3./4.)*its.GetNActiveITSLayers()));
+  printf("Min number of hits requested: %d\n",its.GetMinITSHits());
+  //
+  // max cluster-track chi2 
+  its.RequireMaxChi2Cl(16.);
+  // max chi2/NDF of the Kalman fit
+  its.RequireMaxNormChi2NDF(5.);
+  // penalty to chi2 from missing hit
+  KMCProbe::SetMissingHitPenalty(6);
+  //
+  // Detector, tracking initialization <<<
+
+  its.InitMCWatchHistos();
+  summary.Clear();
+  //
+  int npt = sizeof(pts)/sizeof(double);
+  double ptMin = pts[0];
+  double ptMax = pts[npt-1];
+  //
+  for (int ip=npt;ip--;) {
+    //
+    int nevR = nev;
+    double pt = pts[ip];//1./(ptminI+dpt*ip);  // ptmin+ip*(ptmax-ptmin)/npt;
+
+    double chi2Cl = 16. - (16.-9.)*(1./pt-1./ptMax)/(1./ptMin-1./ptMax);
+    //    if (pt<0.3) nevR = int(0.25*nevR);
+    //    else if (pt<0.5) nevR = int(0.7*nevR);
+    printf("Doing pt(%d)=%.3f for mass %.3f, %d ev | chi2Cl=%.2f\n",ip,pt,mass,nevR,chi2Cl);
+    TString smn = Form("pt%d",ip);
+    TObjArray* trSum = CreateTrackConditions(smn.Data());
+    if (trSum) summary.AddLast(trSum);
+    //
+    its.SolveSingleTrack(mass, pt, eta, trSum, nevR);
+  }
+  //
+  TString smout = sumOut;
+  if (!smout.IsNull()) {
+    gSystem->ExpandPathName(smout);
+    TFile* fl = TFile::Open(smout.Data(),"recreate");
+    if (!fl) {printf("Failed to open %s file\n",smout.Data()); return;}
+    fl->WriteObject(&summary,"trSum","kSingleKey");
+    fl->Close();
+    delete fl;
+  }
+}
+
+//_____________________________________________________________
+TObjArray* CreateTrackConditions(const char* nm)
+{
+  // create set of track criteria to extract the summary
+  KMCTrackSummary *sm = 0;
+  TObjArray* arr = new TObjArray();
+  //
+  int nlr = KMCProbe::GetNITSLayers();
+  //-------------------- put here user code -------------------
+  // summary for ideal correct tracks
+  sm = new KMCTrackSummary("corrAll","corrAll",nlr);
+  sm->SetMinMaxClITS(nlr);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  arr->AddLast(sm);
+  //
+  // summary for good correct tracks
+  sm = new KMCTrackSummary("corr4h","corr4h",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  arr->AddLast(sm);
+  //
+  // summary for good correct tracks with enough points at inner single layers
+  sm = new KMCTrackSummary("corr4hs2in","corr4hs2in",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 layers must have correct cluster
+  sm->AddPatternITSCorr( sm->Bits(0,1));   
+  arr->AddLast(sm);
+  //
+  // summary for good correct tracks with enough points at inner double layers
+  sm = new KMCTrackSummary("corr4hd2in","corr5hd2in",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 double layers must have correct cluster
+  sm->AddPatternITSCorr( sm->Bits(0,1,1));   
+  arr->AddLast(sm);
+  //
+  //
+  // summary for golden tracks with AllNew RS setup: good pattern for offset and pt measurement
+  sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  sm->AddPatternITSCorr( sm->Bits(1,1));    // the innermost 2 double layers must have correct cluster
+  sm->AddPatternITSCorr( sm->Bits(0,0,1)); 
+  //
+  sm->AddPatternITSCorr( sm->Bits(0,0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
+  //
+  sm->AddPatternITSCorr( sm->Bits(0,0,0, 0,0, 1,1)); // measurement at large radius
+  //
+  arr->AddLast(sm);
+  //
+  // summary for golden tracks with curr setup: good pattern for offset and pt measurement
+  sm = new KMCTrackSummary("corr4hdGoodPatt","corr4hdGoodPatt",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(0,0);
+  sm->SetNamePrefix(nm);
+  sm->AddPatternITSCorr( sm->Bits(1,0));    // the innermost 2 double layers must have correct cluster
+  sm->AddPatternITSCorr( sm->Bits(0,1)); 
+  //
+  sm->AddPatternITSCorr( sm->Bits(0,0,1,1)); // measurement in the middle (would be prefereable 2 points!)
+  //
+  sm->AddPatternITSCorr( sm->Bits(0,0, 0,0, 1,1)); // measurement at large radius
+  //
+  arr->AddLast(sm);
+  //
+  // summary for tracks with at least 1 fake cluster
+  sm = new KMCTrackSummary("fakeAny","fakeAny",nlr);
+  sm->SetMinMaxClITS(4);
+  sm->SetMinMaxClITSFake(1);
+  sm->SetNamePrefix(nm);
+  arr->AddLast(sm);
+  //  
+  // // summary for tracks with at least 1 fake cluster (but still acceptable?)
+  // sm = new KMCTrackSummary("corrF","corrF",nlr);
+  // sm->SetMinMaxClITS(4);
+  // sm->SetMinMaxClITSCorr(3);              // at least 3 correct clusters
+  // sm->SetMinMaxClITSFake(1);              // at least 1 fake clusters
+  // //  /*
+  // sm->AddPatternITSCorr( sm->Bits(1));    // the innermost 2 layers must have correct cluster
+  // sm->AddPatternITSCorr( sm->Bits(0,1,1));   
+  // //  */
+  // //  sm->AddPatternITSCorr( sm->Bits(1,1));    // 
+  // // sm->AddPatternITSCorr( sm->Bits(0,0,1,1,1,1));   
+
+  // //  sm->AddPatternITSFakeExcl( sm->Bits(1,1,1)); // no fakes in innermost 3 layers
+  // sm->SetNamePrefix(nm);
+  // arr->AddLast(sm);
+  //  
+  //-----------------------------------------------------------
+  return arr->GetEntries()>0 ? arr : 0;
+}