macro to test ITS TPC matching with FT
authorshahoian <ruben.shahoyan@cern.ch>
Fri, 23 May 2014 18:02:53 +0000 (20:02 +0200)
committershahoian <ruben.shahoyan@cern.ch>
Fri, 23 May 2014 18:02:53 +0000 (20:02 +0200)
ITS/UPGRADE/testITSU/testDetRS.C [new file with mode: 0644]

diff --git a/ITS/UPGRADE/testITSU/testDetRS.C b/ITS/UPGRADE/testITSU/testDetRS.C
new file mode 100644 (file)
index 0000000..c5a905e
--- /dev/null
@@ -0,0 +1,488 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TMath.h>
+#include <TRandom.h>
+#include <TObjArray.h>
+#include "DetectorK.h"
+#include "AliExternalTrackParam.h"
+#endif
+
+const double PTMIN = 0.03;
+const double PTMAX = 30.0;
+const Bool_t VERBOSE = kFALSE;
+const double kMaxChi2 = 10;
+const int    kNbChi2 = 100;
+const int    kNDF = 5;
+
+const double kFactdNdEta = 2;
+
+const UInt_t kFailedBit = 0x1<<14;
+// You have to load the class before ... ;-)
+// .L DetectorK.cxx++
+
+
+const double restR[] = {85.20, 100.20, 120.45,   136.10, 160.10, 220.20, 290};
+
+TrackSol* ts=0;
+DetectorK* its=0;
+
+TH2F* chiComp=0;
+
+TH1F* solPtHisto=0;
+TH1F* effITSSAHisto = 0;
+TH1F* hChiSig = 0;
+TObjArray  arrMatchHisto;
+TObjArray  arrFakeHisto;
+TObjArray  chiBgMatch;
+TObjArray* sigPoolITS=0,*sigPoolGlo=0;
+TObjArray* bgPoolITS=0;
+Double_t wghEtaRange=1;
+
+TObjArray* GenerateSigPool(double ptMin=0.1, double ptMax=20, int nPt=100);
+TObjArray* GenerateBgPool(int nFactor, double ptMin, double ptMax, double refR, double maxDZ, TObjArray* sigPool);
+Double_t   GetRandomPt(double ptMin,double ptMax);
+TH1F*      TestMatching(TrackSol* trGlo, double rLastUpd, double refR, int nPhiTest);
+TObjArray* TestMatchingN(TrackSol* trGlo, const double* rLastUpd, int nRLastUpd, double refR, int nPhiTest);
+TH1F*      PrepChi2Ndf(int ndf=kNDF, int ngen=1000000);
+Double_t   CalcCorrProb(TH1* chiSig, TH1* chiBg, double &smBgRet, double maxChi=-1);
+Bool_t     PassFastCheck(AliExternalTrackParam &tr1, AliExternalTrackParam &tr2, double cut=10);
+
+void testDetRS(double ptMin=0.1, double ptMax=20, int nPt=100, int nFactor=200, double maxDZ=5, double refR=45) 
+{
+  //
+  if (ptMin<PTMIN || ptMin>PTMAX || ptMax<PTMIN || ptMax>PTMAX || ptMin>=ptMax) {
+    printf("Pt range should be within %f %f\n",PTMIN,PTMAX);
+    return;
+  }
+  its = new DetectorK("ALICE","ITS");
+  its->SetAvgRapidity(0.2);
+  //
+  its->AddLayer((char*)"vertex",   0,      0); // dummy vertex for matrix calculation
+  its->AddLayer((char*)"bpipe",  2.0, 0.0022);
+  // new ideal Pixel properties?
+  Double_t x0IB     = 0.003;
+  Double_t x0OB     = 0.008;
+  Double_t resRPhiIB     = 0.0004;
+  Double_t resZIB        = 0.0004;
+  Double_t resRPhiOB     = 0.0004;
+  Double_t resZOB        = 0.0004;
+  Double_t eff           = 0.95;
+  //
+  its->AddLayer((char*)"ddd1",  2.32 ,  x0IB, resRPhiIB, resZIB,eff); 
+  its->AddLayer((char*)"ddd2",  3.13 ,  x0IB, resRPhiIB, resZIB,eff); 
+  its->AddLayer((char*)"ddd3",  3.91 ,  x0IB, resRPhiIB, resZIB,eff); 
+  its->AddLayer((char*)"ddd4",  19.41,  x0OB, resRPhiOB, resZOB,eff); 
+  its->AddLayer((char*)"ddd5",  24.71 ,  x0OB, resRPhiOB, resZOB,eff); 
+  its->AddLayer((char*)"ddd6",  35.33 ,  x0OB, resRPhiOB, resZOB,eff); 
+  its->AddLayer((char*)"ddd7",  40.53 ,  x0OB, resRPhiOB, resZOB,eff); 
+  its->AddLayer((char*)"shell",   43.00 ,  0.01);
+  //
+  printf("Reference radius: %f\n",refR);
+  //
+  its->SetAtLeastHits(5);
+  its->SetAtLeastCorr(5);
+  its->SetAtLeastFake(1);
+  //
+  its->PrintLayout();
+  //
+  sigPoolITS = GenerateSigPool(ptMin,ptMax,nPt);
+  bgPoolITS  = GenerateBgPool(nFactor, ptMin,ptMax, refR, maxDZ, sigPoolITS);
+  //
+  its->AddTPC(0.1,0.1);
+
+  its->AddTRD();
+
+  //
+  sigPoolGlo = GenerateSigPool(ptMin,ptMax,nPt);
+  //
+  //  its->SolveTrack(*ts);
+  //  its->CalcITSEff(*ts,kTRUE);
+  //
+  effITSSAHisto = (TH1F*)solPtHisto->Clone("itsSAeff");
+  effITSSAHisto->SetTitle("ITS SA eff");
+  effITSSAHisto->Reset();
+  hChiSig = PrepChi2Ndf();
+  int nrTest = sizeof(restR)/sizeof(double);
+  for (int ir=0;ir<nrTest;ir++) {
+    TH1F* hmt = (TH1F*)solPtHisto->Clone(Form("its_tpc_match_rupd%.1f_rref%.1f",restR[ir],refR));
+    hmt->SetTitle(Form("ITS-TPC match eff, R TPCmin = %.1f, Rmatch = f%.1f",restR[ir],refR));
+    hmt->Reset();
+    arrMatchHisto.AddLast(hmt);
+    TH1F* hfm = (TH1F*)solPtHisto->Clone(Form("its_tpc_fake_rupd%.1f_rref%.1f",restR[ir],refR));
+    hfm->SetTitle(Form("ITS-TPC fake match eff, R TPCmin = %.1f, Rmatch = f%.1f",restR[ir],refR));
+    hfm->Reset();
+    arrFakeHisto.AddLast(hfm);
+  }
+  //
+  for (int ipt=0;ipt<nPt;ipt++) {
+    ts = (TrackSol*)sigPoolITS->At(ipt);
+    if (ts && !ts->TestBit(kFailedBit)) effITSSAHisto->SetBinContent(ipt+1, ts->fProb[2][0]); // combined SA eff
+    //
+    ts = (TrackSol*)sigPoolGlo->At(ipt);
+    if (ts->TestBit(kFailedBit)) continue;
+    TObjArray* harr = TestMatchingN(ts,restR,nrTest,refR,50);
+    printf("ipt = %d, pt=%.3f : %d histos vs %d tested\n",ipt,ts->fPt, harr ? harr->GetEntries() : 0, nrTest);
+    if (!harr) continue;
+    chiBgMatch.AddAtAndExpand(harr,ipt);
+    //
+    for (int ir=0;ir<nrTest;ir++) {
+      TH1* hbg = (TH1*)harr->At(ir);
+      if (!hbg) continue;
+      double bg = 0;
+      double vl = CalcCorrProb(hChiSig,hbg,bg); // matching prob
+      double effITSSA = effITSSAHisto->GetBinContent(ipt+1);
+      ((TH1*)arrMatchHisto.At(ir))->SetBinContent(ipt+1, vl*effITSSA);
+      ((TH1*)arrFakeHisto.At(ir))->SetBinContent(ipt+1, (1.-vl)*effITSSA + bg*(1.-effITSSA));
+    }
+  }
+  //
+}
+
+TObjArray* GenerateSigPool(double ptMin, double ptMax, int nPt)
+{
+  // evaluate signal tracks
+  TObjArray* arr = new TObjArray(nPt);
+  //
+  if (!solPtHisto) {
+    double dx = log(ptMax/ptMin)/nPt;
+    double *xax = new Double_t[nPt+1];
+    for (int i=0;i<=nPt;i++) xax[i]= ptMin*exp(dx*i);
+    solPtHisto = new TH1F("solPtHisto","solved pt",nPt, xax);
+    delete[] xax;
+  }
+  //
+  for (int ip=0;ip<nPt;ip++) {
+    double pt = solPtHisto->GetBinCenter(ip+1);
+    ts = new TrackSol(its->GetNumberOfLayers(), pt, its->GetAvgRapidity(), -1);
+    if (!its->SolveTrack(*ts) || !its->CalcITSEff(*ts,VERBOSE)) ts->SetBit(kFailedBit);
+    arr->AddAtAndExpand(ts,ip);
+  }
+  //
+  return arr;
+}
+
+TObjArray* GenerateBgPool(int nFactor, double ptMin, double ptMax, double refR, double maxDZ, TObjArray* sigPool)
+{
+  // Generate bg ITS tracks, which may fall within maxDZ of the average 
+  // rapidity track extrapolation at refR
+  //
+  enum {kY,kZ,kSnp,kTgl,kPtI};              // track parameter aliases
+  enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2,kNCovPar}; // cov.matrix aliases
+  //
+  // estimate eta range of generation
+  double thtAv = 2*TMath::ATan(TMath::Exp(-its->GetAvgRapidity()));
+  double zAv = refR*TMath::Tan(thtAv);
+  double zMin = zAv - maxDZ;
+  double zMax = zAv + maxDZ;
+  double thtMin = TMath::ATan(zMin/refR);
+  double thtMax = TMath::ATan(zMax/refR);
+  double etaMax =-TMath::Log( TMath::Tan(thtMin/2.));
+  double etaMin =-TMath::Log( TMath::Tan(thtMax/2.));
+  //
+  wghEtaRange = (etaMax-etaMin); // weight wrt to dn/deta due to the restricted eta range
+  AliExternalTrackParam probTr;
+  double bGauss = its->GetBField()*10;
+  //
+  int ndNdEta = its->GetdNdEtaCent()*kFactdNdEta;
+  int nGen = ndNdEta*nFactor;
+  //
+  int lrRef = its->FindLayerID(refR,0); // compare at this layer
+  //
+  TObjArray* genArr = new TObjArray();
+  //
+  for (int it=0;it<nGen;it++) {
+    double pt = GetRandomPt(ptMin,ptMax);
+    //
+    int ptBin = solPtHisto->FindBin(pt)-1;
+    TrackSol* trSol = (TrackSol*)sigPool->At(ptBin); // fetch the solved track for this pt
+    //
+    if (trSol->TestBit(kFailedBit)) continue;
+    //
+    double eta = etaMin+(etaMax-etaMin)*gRandom->Rndm();
+    int q = gRandom->Rndm()>0.5 ? -1:1;
+    double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta));
+    probTr.Reset();
+    double *trPars = (double*)probTr.GetParameter();
+    double *trCov  = (double*)probTr.GetCovariance();
+    trPars[kY] = 0;                         // start from Y = 0
+    trPars[kZ] = 0;                         //            Z = 0 
+    trPars[kSnp] = 0;                       //            track along X axis at the vertex
+    trPars[kTgl] = TMath::Tan(lambda);      //            dip
+    trPars[kPtI] = q/pt;                    //            q/pt      
+    trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;
+    if (!its->PropagateToR(&probTr,refR,bGauss,1)) {
+      printf("FAILED to propagate to r=%f\n",refR);
+      probTr.Print();
+      continue;
+    }
+    double pos[3];
+    probTr.GetXYZ(pos);  // lab position
+    double phi = TMath::ATan2(pos[1],pos[0]);
+    if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
+    if (!probTr.Rotate(phi))  {
+      printf("FAILED to rotate to phi=%f\n",phi);
+      probTr.Print();
+      continue;
+    }
+    //
+    // assign the error matrix of the ITS outward track with closest pt, estimated with FT
+    AliExternalTrackParam* tr = (AliExternalTrackParam*)trSol->fTrackOutA[lrRef]; // and track pars at requested layer
+    //printf("Pick at lr %d   %p\n",lrRef,tr);
+    //if (tr) tr->Print();
+    memcpy(trCov, tr->GetCovariance(), kNCovPar*sizeof(double));
+    //
+    genArr->AddLast(new AliExternalTrackParam(probTr));
+    solPtHisto->Fill(pt);
+  }
+  //
+  return genArr;
+}
+
+Double_t GetRandomPt(double ptMin,double ptMax)
+{
+  // generate random thermal pt
+  static TF1* dndpt = 0;
+  if (!dndpt) {
+    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);
+    dndpt->SetNpx(1000);
+  }
+  return dndpt->GetRandom(ptMin,ptMax);
+  //
+}
+
+TH1F* TestMatching(TrackSol* trGlo, double rLastUpd, double refR, int nPhiTest)
+{
+  //
+  int maxLrAcc = trGlo->fTrackInw.GetEntries()-1-3; // need at least 3 points fitted
+  Int_t lrLastUpd = its->FindLayerID(rLastUpd, -1); // find the layer below last update R
+  if (lrLastUpd<0 || trGlo->TestBit(kFailedBit)) return 0;
+  //
+  if (lrLastUpd >= maxLrAcc) return 0;
+  // pick the extrapolation to the layer below last requested update layer
+  AliExternalTrackParam* trSignal = (AliExternalTrackParam*)trGlo->fTrackInw[lrLastUpd];
+  if (!trSignal) return 0; 
+  AliExternalTrackParam trSignalC = *trSignal; // use its copy
+  //
+  if (!its->ExtrapolateToR(&trSignalC, refR, trGlo->fMass)) return 0;
+  // bring to the tracking frame of the matching layer
+  //
+  double pos[3];
+  trSignalC.GetXYZ(pos);  // lab position
+  double phi = TMath::ATan2(pos[1],pos[0]);
+  if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
+  if (!trSignalC.Rotate(phi))  {
+    printf("FAILED to rotate to phi=%f\n",phi);
+    trSignalC.Print();
+    return 0;
+  }
+  //
+  int nbg = bgPoolITS->GetEntriesFast();
+  double posBg[3];
+  double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
+  AliExternalTrackParam trBg;
+  //
+  TString tts = Form("hchi_pt%.1f_r%.1f",trSignalC.Pt()*1000,rLastUpd);
+  TH1F* hchi = new TH1F(tts.Data(),tts.Data(),kNbChi2,0,kMaxChi2);
+  hchi->Sumw2();
+  double bGauss = its->GetBField()*10;
+  //
+  for (int ib=0;ib<nbg;ib++) {
+    trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
+    trBg.GetXYZ(posBg);  // lab position
+    //
+    xBg     = trBg.GetX();
+    alpBg   = trBg.GetAlpha();
+    parBg = (double*)trBg.GetParameter();
+    covBg = (double*)trBg.GetCovariance();
+    //
+    memcpy(&parBgOr,parBg,5*sizeof(double));
+    memcpy(&covBgOr,covBg,15*sizeof(double));
+    //
+    // simulate random phi
+    for (int iphi=0;iphi<nPhiTest;iphi++) {
+      double alpBgR = (gRandom->Rndm() - 0.5)*TMath::Pi();
+      trBg.Set(xBg,alpBgR,parBg,covBg);
+      //
+      if (trBg.Rotate(trSignalC.GetAlpha()) && trBg.PropagateTo( trSignalC.GetX(), bGauss)) {
+       if (!PassFastCheck(trSignalC,trBg)) continue;
+       double chi2 = trSignalC.GetPredictedChi2(&trBg);
+       //      trSignalC.Print();
+       //      trBg.Print();
+       //      printf("chi2: %f\n\n\n",chi2);
+       hchi->Fill(chi2/kNDF);
+      }
+      //
+      trBg.Set(xBg,alpBg,parBgOr, covBgOr);  // restore original track
+    }
+    //
+  }
+  //
+  double nTest = nbg*nPhiTest;
+  if (nTest<1) return 0;
+  double nCorr = wghEtaRange*its->GetdNdEtaCent()*kFactdNdEta;
+  printf("Expected: %.0f Tested: %.0f\n",nCorr, nTest);
+  hchi->Scale(nCorr/nTest);
+  return hchi;
+}
+
+//______________________________________________________________________________
+TObjArray* TestMatchingN(TrackSol* trGlo, const double* rLastUpd, int nRLastUpd, double refR, int nPhiTest)
+{
+  //
+  if (trGlo->TestBit(kFailedBit)) return 0;
+  AliExternalTrackParam trSignalC[nRLastUpd];
+  //
+  TObjArray *arrH = new TObjArray(nRLastUpd);
+  int maxLrAcc = trGlo->fTrackInw.GetEntries()-1-3; // need at least 3 points fitted
+  //
+  int sgRef = -1;
+  for (int ir=0;ir<nRLastUpd;ir++) {
+    trSignalC[ir].SetBit(kFailedBit);
+    double rl = rLastUpd[ir];
+    Int_t lrLastUpd = its->FindLayerID(rl, -1); // find the layer below last update R
+    if (lrLastUpd<0) continue;
+    if (lrLastUpd >= maxLrAcc) continue;
+    //
+    // pick the extrapolation to the layer below last requested update layer
+    AliExternalTrackParam* trSignal = (AliExternalTrackParam*)trGlo->fTrackInw[lrLastUpd];
+    if (!trSignal) continue; 
+    trSignalC[ir] = *trSignal; // use its copy
+    //
+    //    printf("\n\n\nPicked at lr=%d r=%f  ",lrLastUpd, rl); trSignal->Print();
+    if (!its->ExtrapolateToR(&trSignalC[ir], refR, trGlo->fMass)) return 0;
+    // bring to the tracking frame of the matching layer
+    //
+    double pos[3];
+    trSignalC[ir].GetXYZ(pos);  // lab position
+    double phi = TMath::ATan2(pos[1],pos[0]);
+    if (!trSignalC[ir].Rotate(phi))  {
+      printf("FAILED to rotate to phi=%f\n",phi);
+      trSignalC[ir].Print();
+      continue;
+    }
+    trSignalC[ir].ResetBit(kFailedBit); // validate the track
+    //
+    //    printf("Taken to r=%f  ",refR); trSignalC[ir].Print();
+    //
+    TString tts = Form("hchi_pt%.1f_r%.1f",trGlo->fPt*1000,rl);
+    TH1F* hchi = new TH1F(tts.Data(),tts.Data(),kNbChi2,0,kMaxChi2);
+    hchi->Sumw2();
+    arrH->AddAt(hchi,ir);
+    if (sgRef<0) sgRef = ir;
+  }
+  //
+  if (sgRef<0) return 0;
+  //
+  int nbg = bgPoolITS->GetEntriesFast();
+  double posBg[3];
+  double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
+  AliExternalTrackParam trBg;
+  //
+  double bGauss = its->GetBField()*10;
+  //
+  for (int ib=0;ib<nbg;ib++) {
+    trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
+    trBg.GetXYZ(posBg);  // lab position
+    //
+    xBg     = trBg.GetX();
+    alpBg   = trBg.GetAlpha();
+    parBg = (double*)trBg.GetParameter();
+    covBg = (double*)trBg.GetCovariance();
+    //
+    memcpy(&parBgOr,parBg,5*sizeof(double));
+    memcpy(&covBgOr,covBg,15*sizeof(double));
+    //
+    // simulate random phi
+    for (int iphi=0;iphi<nPhiTest;iphi++) {
+      double alpBgR = (gRandom->Rndm() - 0.5)*TMath::Pi();
+      trBg.Set(xBg,alpBgR,parBg,covBg);
+      //
+      if (trBg.Rotate(trSignalC[sgRef].GetAlpha()) && trBg.PropagateTo( trSignalC[sgRef].GetX(), bGauss)) {
+       for (int ir=0;ir<nRLastUpd;ir++) {
+         if (trSignalC[ir].TestBit(kFailedBit)) continue;
+         if (!PassFastCheck(trSignalC[ir],trBg)) continue;
+         double chi2 = trSignalC[ir].GetPredictedChi2(&trBg);
+         ((TH1*)arrH->UncheckedAt(ir))->Fill(chi2/kNDF);
+       }
+       //
+      }
+      trBg.Set(xBg,alpBg,parBgOr, covBgOr);  // restore original track
+      //
+    }
+  }
+  //
+  double nTest = nbg*nPhiTest;
+  if (nTest<1) return 0;
+  double nCorr = wghEtaRange*its->GetdNdEtaCent()*kFactdNdEta;
+  printf("Expected: %.0f Tested: %.0f\n",nCorr, nTest);
+  for (int ir=0;ir<nRLastUpd;ir++) {
+    TH1* hchi = (TH1*)arrH->UncheckedAt(ir);
+    if (!hchi) continue;
+    hchi->Scale(nCorr/nTest);
+  }
+  return arrH;
+}
+
+//_________________________________________________________________________________________
+TH1F* PrepChi2Ndf(int ndf, int ngen)
+{
+  // generate distribution for chi2 with given ndf
+  TString ttl = Form("Chi2_NDF%d",ndf);
+  TH1F* hh = new TH1F(ttl.Data(),ttl.Data(),kNbChi2,0,kMaxChi2);
+  for (int i=0;i<ngen;i++) hh->Fill(TMath::ChisquareQuantile(gRandom->Rndm(),ndf)/ndf);
+  hh->Scale(1./ngen);
+  return hh;
+}
+
+//____________________________________________________
+Double_t CalcCorrProb(TH1* chiSig, TH1* chiBg, double &smBgRet, double maxChi) 
+{
+  // Calculate probability of all fake matches being worse than correct one
+  // The chi2 histos should be normalized
+  int nb = chiSig->GetNbinsX();
+  if (maxChi>0) {
+    int nt = chiSig->FindBin(maxChi);
+    if (nt<nb) nb = nt;
+  }
+  //
+  double bgchiC = maxChi;
+  if (bgchiC<0) bgchiC = 5.;
+  int maxBgBinCut = chiSig->FindBin(bgchiC);
+  double smBg = 0;
+  smBgRet = 0;
+  double probTotCorr = 0;
+  for (int i=1;i<=nb;i++) {
+    double wSig = chiSig->GetBinContent(i);
+    probTotCorr += wSig*TMath::Exp(-smBg);  // Poisson prob of not having a bg with better chi2
+    smBg += chiBg->GetBinContent(i);
+    if (i<=maxBgBinCut) smBgRet = smBg;
+  }
+  //
+  smBg = 1.-TMath::Exp(smBg); // prob. to find bg within a cut in absence of correct match
+  return probTotCorr;
+}
+
+//__________________________________________________
+Bool_t PassFastCheck(AliExternalTrackParam &tr1, AliExternalTrackParam &tr2, double cut)
+{
+  enum {kY,kZ,kSnp,kTgl,kPtI};              // track parameter aliases
+  enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
+  const int erId[5] = {kY2,kZ2,kSnp2,kTgl2,kPtI2};
+  //
+  // fast check ignoring non-diagonal elements
+  double *cov1 = (double*)tr1.GetCovariance();
+  double *cov2 = (double*)tr2.GetCovariance();
+  double *par1 = (double*)tr1.GetParameter();
+  double *par2 = (double*)tr2.GetParameter();
+  //
+  double chi2 = 0;
+  for (int i=0;i<5;i++) {
+    double df = par1[i]-par2[i];
+    double err2 = cov1[erId[i]]+cov2[erId[i]];
+    chi2 += df*df/err2;
+  }
+  return (chi2>cut*kNDF) ? kFALSE : kTRUE;
+  //
+}