ITS/UPGRADE/KMCDetector.cxx

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