ITS/UPGRADE/testITSU/testDetRS.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TH1.h>
   3 #include <TH2.h>
   4 #include <TF1.h>
   5 #include <TMath.h>
   6 #include <TRandom.h>
   7 #include <TObjArray.h>
   8 #include "DetectorK.h"
   9 #include "AliExternalTrackParam.h"
  10 #endif
  11
  12 const double PTMIN = 0.03;
  13 const double PTMAX = 30.0;
  14 const Bool_t VERBOSE = kFALSE;
  15 const double kMaxChi2 = 10;
  16 const int    kNbChi2 = 100;
  17 const int    kNDF = 5;
  18
  19 const double kFactdNdEta = 2;
  20
  21 const UInt_t kFailedBit = 0x1<<14;
  22 // You have to load the class before ... ;-)
  23 // .L DetectorK.cxx++
  24
  25
  26 const double restR[] = {85.20, 100.20, 120.45,   136.10, 160.10, 220.20, 290};
  27
  28 TrackSol* ts=0;
  29 DetectorK* its=0;
  30
  31 TH2F* chiComp=0;
  32
  33 TH1F* solPtHisto=0;
  34 TH1F* effITSSAHisto = 0;
  35 TH1F* hChiSig = 0;
  36 TObjArray  arrMatchHisto;
  37 TObjArray  arrFakeHisto;
  38 TObjArray  chiBgMatch;
  39 TObjArray* sigPoolITS=0,*sigPoolGlo=0;
  40 TObjArray* bgPoolITS=0;
  41 Double_t wghEtaRange=1;
  42
  43 TObjArray* GenerateSigPool(double ptMin=0.1, double ptMax=20, int nPt=100);
  44 TObjArray* GenerateBgPool(int nFactor, double ptMin, double ptMax, double refR, double maxDZ, TObjArray* sigPool);
  45 Double_t   GetRandomPt(double ptMin,double ptMax);
  46 TH1F*      TestMatching(TrackSol* trGlo, double rLastUpd, double refR, int nPhiTest);
  47 TObjArray* TestMatchingN(TrackSol* trGlo, const double* rLastUpd, int nRLastUpd, double refR, int nPhiTest);
  48 TH1F*      PrepChi2Ndf(int ndf=kNDF, int ngen=1000000);
  49 Double_t   CalcCorrProb(TH1* chiSig, TH1* chiBg, double &smBgRet, double maxChi=-1);
  50 Bool_t     PassFastCheck(AliExternalTrackParam &tr1, AliExternalTrackParam &tr2, double cut=10);
  51
  52 void testDetRS(double ptMin=0.1, double ptMax=20, int nPt=100, int nFactor=200, double maxDZ=5, double refR=45)
  53 {
  54   //
  55   if (ptMin<PTMIN || ptMin>PTMAX || ptMax<PTMIN || ptMax>PTMAX || ptMin>=ptMax) {
  56     printf("Pt range should be within %f %f\n",PTMIN,PTMAX);
  57     return;
  58   }
  59   its = new DetectorK("ALICE","ITS");
  60   its->SetAvgRapidity(0.2);
  61   //
  62   its->AddLayer((char*)"vertex",   0,      0); // dummy vertex for matrix calculation
  63   its->AddLayer((char*)"bpipe",  2.0, 0.0022);
  64   // new ideal Pixel properties?
  65   Double_t x0IB     = 0.003;
  66   Double_t x0OB     = 0.008;
  67   Double_t resRPhiIB     = 0.0004;
  68   Double_t resZIB        = 0.0004;
  69   Double_t resRPhiOB     = 0.0004;
  70   Double_t resZOB        = 0.0004;
  71   Double_t eff           = 0.95;
  72   //
  73   its->AddLayer((char*)"ddd1",  2.32 ,  x0IB, resRPhiIB, resZIB,eff);
  74   its->AddLayer((char*)"ddd2",  3.13 ,  x0IB, resRPhiIB, resZIB,eff);
  75   its->AddLayer((char*)"ddd3",  3.91 ,  x0IB, resRPhiIB, resZIB,eff);
  76   its->AddLayer((char*)"ddd4",  19.41,  x0OB, resRPhiOB, resZOB,eff);
  77   its->AddLayer((char*)"ddd5",  24.71 ,  x0OB, resRPhiOB, resZOB,eff);
  78   its->AddLayer((char*)"ddd6",  35.33 ,  x0OB, resRPhiOB, resZOB,eff);
  79   its->AddLayer((char*)"ddd7",  40.53 ,  x0OB, resRPhiOB, resZOB,eff);
  80   its->AddLayer((char*)"shell",   43.00 ,  0.01);
  81   //
  82   printf("Reference radius: %f\n",refR);
  83   //
  84   its->SetAtLeastHits(5);
  85   its->SetAtLeastCorr(5);
  86   its->SetAtLeastFake(1);
  87   //
  88   its->PrintLayout();
  89   //
  90   sigPoolITS = GenerateSigPool(ptMin,ptMax,nPt);
  91   bgPoolITS  = GenerateBgPool(nFactor, ptMin,ptMax, refR, maxDZ, sigPoolITS);
  92   //
  93   its->AddTPC(0.1,0.1);
  94
  95   its->AddTRD();
  96
  97   //
  98   sigPoolGlo = GenerateSigPool(ptMin,ptMax,nPt);
  99   //
 100   //  its->SolveTrack(*ts);
 101   //  its->CalcITSEff(*ts,kTRUE);
 102   //
 103   effITSSAHisto = (TH1F*)solPtHisto->Clone("itsSAeff");
 104   effITSSAHisto->SetTitle("ITS SA eff");
 105   effITSSAHisto->Reset();
 106   hChiSig = PrepChi2Ndf();
 107   int nrTest = sizeof(restR)/sizeof(double);
 108   for (int ir=0;ir<nrTest;ir++) {
 109     TH1F* hmt = (TH1F*)solPtHisto->Clone(Form("its_tpc_match_rupd%.1f_rref%.1f",restR[ir],refR));
 110     hmt->SetTitle(Form("ITS-TPC match eff, R TPCmin = %.1f, Rmatch = f%.1f",restR[ir],refR));
 111     hmt->Reset();
 112     arrMatchHisto.AddLast(hmt);
 113     TH1F* hfm = (TH1F*)solPtHisto->Clone(Form("its_tpc_fake_rupd%.1f_rref%.1f",restR[ir],refR));
 114     hfm->SetTitle(Form("ITS-TPC fake match eff, R TPCmin = %.1f, Rmatch = f%.1f",restR[ir],refR));
 115     hfm->Reset();
 116     arrFakeHisto.AddLast(hfm);
 117   }
 118   //
 119   for (int ipt=0;ipt<nPt;ipt++) {
 120     ts = (TrackSol*)sigPoolITS->At(ipt);
 121     if (ts && !ts->TestBit(kFailedBit)) effITSSAHisto->SetBinContent(ipt+1, ts->fProb[2][0]); // combined SA eff
 122     //
 123     ts = (TrackSol*)sigPoolGlo->At(ipt);
 124     if (ts->TestBit(kFailedBit)) continue;
 125     TObjArray* harr = TestMatchingN(ts,restR,nrTest,refR,50);
 126     printf("ipt = %d, pt=%.3f : %d histos vs %d tested\n",ipt,ts->fPt, harr ? harr->GetEntries() : 0, nrTest);
 127     if (!harr) continue;
 128     chiBgMatch.AddAtAndExpand(harr,ipt);
 129     //
 130     for (int ir=0;ir<nrTest;ir++) {
 131       TH1* hbg = (TH1*)harr->At(ir);
 132       if (!hbg) continue;
 133       double bg = 0;
 134       double vl = CalcCorrProb(hChiSig,hbg,bg); // matching prob
 135       double effITSSA = effITSSAHisto->GetBinContent(ipt+1);
 136       ((TH1*)arrMatchHisto.At(ir))->SetBinContent(ipt+1, vl*effITSSA);
 137       ((TH1*)arrFakeHisto.At(ir))->SetBinContent(ipt+1, (1.-vl)*effITSSA + bg*(1.-effITSSA));
 138     }
 139   }
 140   //
 141 }
 142
 143 TObjArray* GenerateSigPool(double ptMin, double ptMax, int nPt)
 144 {
 145   // evaluate signal tracks
 146   TObjArray* arr = new TObjArray(nPt);
 147   //
 148   if (!solPtHisto) {
 149     double dx = log(ptMax/ptMin)/nPt;
 150     double *xax = new Double_t[nPt+1];
 151     for (int i=0;i<=nPt;i++) xax[i]= ptMin*exp(dx*i);
 152     solPtHisto = new TH1F("solPtHisto","solved pt",nPt, xax);
 153     delete[] xax;
 154   }
 155   //
 156   for (int ip=0;ip<nPt;ip++) {
 157     double pt = solPtHisto->GetBinCenter(ip+1);
 158     ts = new TrackSol(its->GetNumberOfLayers(), pt, its->GetAvgRapidity(), -1);
 159     if (!its->SolveTrack(*ts) || !its->CalcITSEff(*ts,VERBOSE)) ts->SetBit(kFailedBit);
 160     arr->AddAtAndExpand(ts,ip);
 161   }
 162   //
 163   return arr;
 164 }
 165
 166 TObjArray* GenerateBgPool(int nFactor, double ptMin, double ptMax, double refR, double maxDZ, TObjArray* sigPool)
 167 {
 168   // Generate bg ITS tracks, which may fall within maxDZ of the average
 169   // rapidity track extrapolation at refR
 170   //
 171   enum {kY,kZ,kSnp,kTgl,kPtI};              // track parameter aliases
 172   enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2,kNCovPar}; // cov.matrix aliases
 173   //
 174   // estimate eta range of generation
 175   double thtAv = 2*TMath::ATan(TMath::Exp(-its->GetAvgRapidity()));
 176   double zAv = refR*TMath::Tan(thtAv);
 177   double zMin = zAv - maxDZ;
 178   double zMax = zAv + maxDZ;
 179   double thtMin = TMath::ATan(zMin/refR);
 180   double thtMax = TMath::ATan(zMax/refR);
 181   double etaMax =-TMath::Log( TMath::Tan(thtMin/2.));
 182   double etaMin =-TMath::Log( TMath::Tan(thtMax/2.));
 183   //
 184   wghEtaRange = (etaMax-etaMin); // weight wrt to dn/deta due to the restricted eta range
 185   AliExternalTrackParam probTr;
 186   double bGauss = its->GetBField()*10;
 187   //
 188   int ndNdEta = its->GetdNdEtaCent()*kFactdNdEta;
 189   int nGen = ndNdEta*nFactor;
 190   //
 191   int lrRef = its->FindLayerID(refR,0); // compare at this layer
 192   //
 193   TObjArray* genArr = new TObjArray();
 194   //
 195   for (int it=0;it<nGen;it++) {
 196     double pt = GetRandomPt(ptMin,ptMax);
 197     //
 198     int ptBin = solPtHisto->FindBin(pt)-1;
 199     TrackSol* trSol = (TrackSol*)sigPool->At(ptBin); // fetch the solved track for this pt
 200     //
 201     if (trSol->TestBit(kFailedBit)) continue;
 202     //
 203     double eta = etaMin+(etaMax-etaMin)*gRandom->Rndm();
 204     int q = gRandom->Rndm()>0.5 ? -1:1;
 205     double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta));
 206     probTr.Reset();
 207     double *trPars = (double*)probTr.GetParameter();
 208     double *trCov  = (double*)probTr.GetCovariance();
 209     trPars[kY] = 0;                         // start from Y = 0
 210     trPars[kZ] = 0;                         //            Z = 0
 211     trPars[kSnp] = 0;                       //            track along X axis at the vertex
 212     trPars[kTgl] = TMath::Tan(lambda);      //            dip
 213     trPars[kPtI] = q/pt;                    //            q/pt
 214     trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;
 215     if (!its->PropagateToR(&probTr,refR,bGauss,1)) {
 216       printf("FAILED to propagate to r=%f\n",refR);
 217       probTr.Print();
 218       continue;
 219     }
 220     double pos[3];
 221     probTr.GetXYZ(pos);  // lab position
 222     double phi = TMath::ATan2(pos[1],pos[0]);
 223     if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
 224     if (!probTr.Rotate(phi))  {
 225       printf("FAILED to rotate to phi=%f\n",phi);
 226       probTr.Print();
 227       continue;
 228     }
 229     //
 230     // assign the error matrix of the ITS outward track with closest pt, estimated with FT
 231     AliExternalTrackParam* tr = (AliExternalTrackParam*)trSol->fTrackOutA[lrRef]; // and track pars at requested layer
 232     //printf("Pick at lr %d   %p\n",lrRef,tr);
 233     //if (tr) tr->Print();
 234     memcpy(trCov, tr->GetCovariance(), kNCovPar*sizeof(double));
 235     //
 236     genArr->AddLast(new AliExternalTrackParam(probTr));
 237     solPtHisto->Fill(pt);
 238   }
 239   //
 240   return genArr;
 241 }
 242
 243 Double_t GetRandomPt(double ptMin,double ptMax)
 244 {
 245   // generate random thermal pt
 246   static TF1* dndpt = 0;
 247   if (!dndpt) {
 248     dndpt = new TF1("dndpt","pow(sqrt(x*x+0.14*0.14),1.5)*exp(-sqrt(x*x+0.14*0.14)/0.17)",PTMIN,PTMAX);
 249     dndpt->SetNpx(1000);
 250   }
 251   return dndpt->GetRandom(ptMin,ptMax);
 252   //
 253 }
 254
 255 TH1F* TestMatching(TrackSol* trGlo, double rLastUpd, double refR, int nPhiTest)
 256 {
 257   //
 258   int maxLrAcc = trGlo->fTrackInw.GetEntries()-1-3; // need at least 3 points fitted
 259   Int_t lrLastUpd = its->FindLayerID(rLastUpd, -1); // find the layer below last update R
 260   if (lrLastUpd<0 || trGlo->TestBit(kFailedBit)) return 0;
 261   //
 262   if (lrLastUpd >= maxLrAcc) return 0;
 263   // pick the extrapolation to the layer below last requested update layer
 264   AliExternalTrackParam* trSignal = (AliExternalTrackParam*)trGlo->fTrackInw[lrLastUpd];
 265   if (!trSignal) return 0;
 266   AliExternalTrackParam trSignalC = *trSignal; // use its copy
 267   //
 268   if (!its->ExtrapolateToR(&trSignalC, refR, trGlo->fMass)) return 0;
 269   // bring to the tracking frame of the matching layer
 270   //
 271   double pos[3];
 272   trSignalC.GetXYZ(pos);  // lab position
 273   double phi = TMath::ATan2(pos[1],pos[0]);
 274   if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
 275   if (!trSignalC.Rotate(phi))  {
 276     printf("FAILED to rotate to phi=%f\n",phi);
 277     trSignalC.Print();
 278     return 0;
 279   }
 280   //
 281   int nbg = bgPoolITS->GetEntriesFast();
 282   double posBg[3];
 283   double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
 284   AliExternalTrackParam trBg;
 285   //
 286   TString tts = Form("hchi_pt%.1f_r%.1f",trSignalC.Pt()*1000,rLastUpd);
 287   TH1F* hchi = new TH1F(tts.Data(),tts.Data(),kNbChi2,0,kMaxChi2);
 288   hchi->Sumw2();
 289   double bGauss = its->GetBField()*10;
 290   //
 291   for (int ib=0;ib<nbg;ib++) {
 292     trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
 293     trBg.GetXYZ(posBg);  // lab position
 294     //
 295     xBg     = trBg.GetX();
 296     alpBg   = trBg.GetAlpha();
 297     parBg = (double*)trBg.GetParameter();
 298     covBg = (double*)trBg.GetCovariance();
 299     //
 300     memcpy(&parBgOr,parBg,5*sizeof(double));
 301     memcpy(&covBgOr,covBg,15*sizeof(double));
 302     //
 303     // simulate random phi
 304     for (int iphi=0;iphi<nPhiTest;iphi++) {
 305       double alpBgR = (gRandom->Rndm() - 0.5)*TMath::Pi();
 306       trBg.Set(xBg,alpBgR,parBg,covBg);
 307       //
 308       if (trBg.Rotate(trSignalC.GetAlpha()) && trBg.PropagateTo( trSignalC.GetX(), bGauss)) {
 309         if (!PassFastCheck(trSignalC,trBg)) continue;
 310         double chi2 = trSignalC.GetPredictedChi2(&trBg);
 311         //      trSignalC.Print();
 312         //      trBg.Print();
 313         //      printf("chi2: %f\n\n\n",chi2);
 314         hchi->Fill(chi2/kNDF);
 315       }
 316       //
 317       trBg.Set(xBg,alpBg,parBgOr, covBgOr);  // restore original track
 318     }
 319     //
 320   }
 321   //
 322   double nTest = nbg*nPhiTest;
 323   if (nTest<1) return 0;
 324   double nCorr = wghEtaRange*its->GetdNdEtaCent()*kFactdNdEta;
 325   printf("Expected: %.0f Tested: %.0f\n",nCorr, nTest);
 326   hchi->Scale(nCorr/nTest);
 327   return hchi;
 328 }
 329
 330 //______________________________________________________________________________
 331 TObjArray* TestMatchingN(TrackSol* trGlo, const double* rLastUpd, int nRLastUpd, double refR, int nPhiTest)
 332 {
 333   //
 334   if (trGlo->TestBit(kFailedBit)) return 0;
 335   AliExternalTrackParam trSignalC[nRLastUpd];
 336   //
 337   TObjArray *arrH = new TObjArray(nRLastUpd);
 338   int maxLrAcc = trGlo->fTrackInw.GetEntries()-1-3; // need at least 3 points fitted
 339   //
 340   int sgRef = -1;
 341   for (int ir=0;ir<nRLastUpd;ir++) {
 342     trSignalC[ir].SetBit(kFailedBit);
 343     double rl = rLastUpd[ir];
 344     Int_t lrLastUpd = its->FindLayerID(rl, -1); // find the layer below last update R
 345     if (lrLastUpd<0) continue;
 346     if (lrLastUpd >= maxLrAcc) continue;
 347     //
 348     // pick the extrapolation to the layer below last requested update layer
 349     AliExternalTrackParam* trSignal = (AliExternalTrackParam*)trGlo->fTrackInw[lrLastUpd];
 350     if (!trSignal) continue;
 351     trSignalC[ir] = *trSignal; // use its copy
 352     //
 353     //    printf("\n\n\nPicked at lr=%d r=%f  ",lrLastUpd, rl); trSignal->Print();
 354     if (!its->ExtrapolateToR(&trSignalC[ir], refR, trGlo->fMass)) return 0;
 355     // bring to the tracking frame of the matching layer
 356     //
 357     double pos[3];
 358     trSignalC[ir].GetXYZ(pos);  // lab position
 359     double phi = TMath::ATan2(pos[1],pos[0]);
 360     if (!trSignalC[ir].Rotate(phi))  {
 361       printf("FAILED to rotate to phi=%f\n",phi);
 362       trSignalC[ir].Print();
 363       continue;
 364     }
 365     trSignalC[ir].ResetBit(kFailedBit); // validate the track
 366     //
 367     //    printf("Taken to r=%f  ",refR); trSignalC[ir].Print();
 368     //
 369     TString tts = Form("hchi_pt%.1f_r%.1f",trGlo->fPt*1000,rl);
 370     TH1F* hchi = new TH1F(tts.Data(),tts.Data(),kNbChi2,0,kMaxChi2);
 371     hchi->Sumw2();
 372     arrH->AddAt(hchi,ir);
 373     if (sgRef<0) sgRef = ir;
 374   }
 375   //
 376   if (sgRef<0) return 0;
 377   //
 378   int nbg = bgPoolITS->GetEntriesFast();
 379   double posBg[3];
 380   double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
 381   AliExternalTrackParam trBg;
 382   //
 383   double bGauss = its->GetBField()*10;
 384   //
 385   for (int ib=0;ib<nbg;ib++) {
 386     trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
 387     trBg.GetXYZ(posBg);  // lab position
 388     //
 389     xBg     = trBg.GetX();
 390     alpBg   = trBg.GetAlpha();
 391     parBg = (double*)trBg.GetParameter();
 392     covBg = (double*)trBg.GetCovariance();
 393     //
 394     memcpy(&parBgOr,parBg,5*sizeof(double));
 395     memcpy(&covBgOr,covBg,15*sizeof(double));
 396     //
 397     // simulate random phi
 398     for (int iphi=0;iphi<nPhiTest;iphi++) {
 399       double alpBgR = (gRandom->Rndm() - 0.5)*TMath::Pi();
 400       trBg.Set(xBg,alpBgR,parBg,covBg);
 401       //
 402       if (trBg.Rotate(trSignalC[sgRef].GetAlpha()) && trBg.PropagateTo( trSignalC[sgRef].GetX(), bGauss)) {
 403         for (int ir=0;ir<nRLastUpd;ir++) {
 404           if (trSignalC[ir].TestBit(kFailedBit)) continue;
 405           if (!PassFastCheck(trSignalC[ir],trBg)) continue;
 406           double chi2 = trSignalC[ir].GetPredictedChi2(&trBg);
 407           ((TH1*)arrH->UncheckedAt(ir))->Fill(chi2/kNDF);
 408         }
 409         //
 410       }
 411       trBg.Set(xBg,alpBg,parBgOr, covBgOr);  // restore original track
 412       //
 413     }
 414   }
 415   //
 416   double nTest = nbg*nPhiTest;
 417   if (nTest<1) return 0;
 418   double nCorr = wghEtaRange*its->GetdNdEtaCent()*kFactdNdEta;
 419   printf("Expected: %.0f Tested: %.0f\n",nCorr, nTest);
 420   for (int ir=0;ir<nRLastUpd;ir++) {
 421     TH1* hchi = (TH1*)arrH->UncheckedAt(ir);
 422     if (!hchi) continue;
 423     hchi->Scale(nCorr/nTest);
 424   }
 425   return arrH;
 426 }
 427
 428 //_________________________________________________________________________________________
 429 TH1F* PrepChi2Ndf(int ndf, int ngen)
 430 {
 431   // generate distribution for chi2 with given ndf
 432   TString ttl = Form("Chi2_NDF%d",ndf);
 433   TH1F* hh = new TH1F(ttl.Data(),ttl.Data(),kNbChi2,0,kMaxChi2);
 434   for (int i=0;i<ngen;i++) hh->Fill(TMath::ChisquareQuantile(gRandom->Rndm(),ndf)/ndf);
 435   hh->Scale(1./ngen);
 436   return hh;
 437 }
 438
 439 //____________________________________________________
 440 Double_t CalcCorrProb(TH1* chiSig, TH1* chiBg, double &smBgRet, double maxChi)
 441 {
 442   // Calculate probability of all fake matches being worse than correct one
 443   // The chi2 histos should be normalized
 444   int nb = chiSig->GetNbinsX();
 445   if (maxChi>0) {
 446     int nt = chiSig->FindBin(maxChi);
 447     if (nt<nb) nb = nt;
 448   }
 449   //
 450   double bgchiC = maxChi;
 451   if (bgchiC<0) bgchiC = 5.;
 452   int maxBgBinCut = chiSig->FindBin(bgchiC);
 453   double smBg = 0;
 454   smBgRet = 0;
 455   double probTotCorr = 0;
 456   for (int i=1;i<=nb;i++) {
 457     double wSig = chiSig->GetBinContent(i);
 458     probTotCorr += wSig*TMath::Exp(-smBg);  // Poisson prob of not having a bg with better chi2
 459     smBg += chiBg->GetBinContent(i);
 460     if (i<=maxBgBinCut) smBgRet = smBg;
 461   }
 462   //
 463   smBg = 1.-TMath::Exp(smBg); // prob. to find bg within a cut in absence of correct match
 464   return probTotCorr;
 465 }
 466
 467 //__________________________________________________
 468 Bool_t PassFastCheck(AliExternalTrackParam &tr1, AliExternalTrackParam &tr2, double cut)
 469 {
 470   enum {kY,kZ,kSnp,kTgl,kPtI};              // track parameter aliases
 471   enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
 472   const int erId[5] = {kY2,kZ2,kSnp2,kTgl2,kPtI2};
 473   //
 474   // fast check ignoring non-diagonal elements
 475   double *cov1 = (double*)tr1.GetCovariance();
 476   double *cov2 = (double*)tr2.GetCovariance();
 477   double *par1 = (double*)tr1.GetParameter();
 478   double *par2 = (double*)tr2.GetParameter();
 479   //
 480   double chi2 = 0;
 481   for (int i=0;i<5;i++) {
 482     double df = par1[i]-par2[i];
 483     double err2 = cov1[erId[i]]+cov2[erId[i]];
 484     chi2 += df*df/err2;
 485   }
 486   return (chi2>cut*kNDF) ? kFALSE : kTRUE;
 487   //
 488 }