1 #if !defined(__CINT__) || defined(__MAKECINT__)
9 #include "AliExternalTrackParam.h"
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;
19 const double kFactdNdEta = 2;
21 const UInt_t kFailedBit = 0x1<<14;
22 // You have to load the class before ... ;-)
26 const double restR[] = {85.20, 100.20, 120.45, 136.10, 160.10, 220.20, 290};
34 TH1F* effITSSAHisto = 0;
36 TObjArray arrMatchHisto;
37 TObjArray arrFakeHisto;
39 TObjArray* sigPoolITS=0,*sigPoolGlo=0;
40 TObjArray* bgPoolITS=0;
41 Double_t wghEtaRange=1;
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);
52 void testDetRS(double ptMin=0.1, double ptMax=20, int nPt=100, int nFactor=200, double maxDZ=5, double refR=45)
55 if (ptMin<PTMIN || ptMin>PTMAX || ptMax<PTMIN || ptMax>PTMAX || ptMin>=ptMax) {
56 printf("Pt range should be within %f %f\n",PTMIN,PTMAX);
59 its = new DetectorK("ALICE","ITS");
60 its->SetAvgRapidity(0.2);
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;
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);
82 printf("Reference radius: %f\n",refR);
84 its->SetAtLeastHits(5);
85 its->SetAtLeastCorr(5);
86 its->SetAtLeastFake(1);
90 sigPoolITS = GenerateSigPool(ptMin,ptMax,nPt);
91 bgPoolITS = GenerateBgPool(nFactor, ptMin,ptMax, refR, maxDZ, sigPoolITS);
98 sigPoolGlo = GenerateSigPool(ptMin,ptMax,nPt);
100 // its->SolveTrack(*ts);
101 // its->CalcITSEff(*ts,kTRUE);
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));
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));
116 arrFakeHisto.AddLast(hfm);
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
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);
128 chiBgMatch.AddAtAndExpand(harr,ipt);
130 for (int ir=0;ir<nrTest;ir++) {
131 TH1* hbg = (TH1*)harr->At(ir);
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));
143 TObjArray* GenerateSigPool(double ptMin, double ptMax, int nPt)
145 // evaluate signal tracks
146 TObjArray* arr = new TObjArray(nPt);
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);
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);
166 TObjArray* GenerateBgPool(int nFactor, double ptMin, double ptMax, double refR, double maxDZ, TObjArray* sigPool)
168 // Generate bg ITS tracks, which may fall within maxDZ of the average
169 // rapidity track extrapolation at refR
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
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.));
184 wghEtaRange = (etaMax-etaMin); // weight wrt to dn/deta due to the restricted eta range
185 AliExternalTrackParam probTr;
186 double bGauss = its->GetBField()*10;
188 int ndNdEta = its->GetdNdEtaCent()*kFactdNdEta;
189 int nGen = ndNdEta*nFactor;
191 int lrRef = its->FindLayerID(refR,0); // compare at this layer
193 TObjArray* genArr = new TObjArray();
195 for (int it=0;it<nGen;it++) {
196 double pt = GetRandomPt(ptMin,ptMax);
198 int ptBin = solPtHisto->FindBin(pt)-1;
199 TrackSol* trSol = (TrackSol*)sigPool->At(ptBin); // fetch the solved track for this pt
201 if (trSol->TestBit(kFailedBit)) continue;
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));
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);
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);
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));
236 genArr->AddLast(new AliExternalTrackParam(probTr));
237 solPtHisto->Fill(pt);
243 Double_t GetRandomPt(double ptMin,double ptMax)
245 // generate random thermal pt
246 static TF1* dndpt = 0;
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);
251 return dndpt->GetRandom(ptMin,ptMax);
255 TH1F* TestMatching(TrackSol* trGlo, double rLastUpd, double refR, int nPhiTest)
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;
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
268 if (!its->ExtrapolateToR(&trSignalC, refR, trGlo->fMass)) return 0;
269 // bring to the tracking frame of the matching layer
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);
281 int nbg = bgPoolITS->GetEntriesFast();
283 double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
284 AliExternalTrackParam trBg;
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);
289 double bGauss = its->GetBField()*10;
291 for (int ib=0;ib<nbg;ib++) {
292 trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
293 trBg.GetXYZ(posBg); // lab position
296 alpBg = trBg.GetAlpha();
297 parBg = (double*)trBg.GetParameter();
298 covBg = (double*)trBg.GetCovariance();
300 memcpy(&parBgOr,parBg,5*sizeof(double));
301 memcpy(&covBgOr,covBg,15*sizeof(double));
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);
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();
313 // printf("chi2: %f\n\n\n",chi2);
314 hchi->Fill(chi2/kNDF);
317 trBg.Set(xBg,alpBg,parBgOr, covBgOr); // restore original track
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);
330 //______________________________________________________________________________
331 TObjArray* TestMatchingN(TrackSol* trGlo, const double* rLastUpd, int nRLastUpd, double refR, int nPhiTest)
334 if (trGlo->TestBit(kFailedBit)) return 0;
335 AliExternalTrackParam trSignalC[nRLastUpd];
337 TObjArray *arrH = new TObjArray(nRLastUpd);
338 int maxLrAcc = trGlo->fTrackInw.GetEntries()-1-3; // need at least 3 points fitted
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;
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
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
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();
365 trSignalC[ir].ResetBit(kFailedBit); // validate the track
367 // printf("Taken to r=%f ",refR); trSignalC[ir].Print();
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);
372 arrH->AddAt(hchi,ir);
373 if (sgRef<0) sgRef = ir;
376 if (sgRef<0) return 0;
378 int nbg = bgPoolITS->GetEntriesFast();
380 double xBg,alpBg,parBgOr[5],covBgOr[15],*parBg,*covBg;
381 AliExternalTrackParam trBg;
383 double bGauss = its->GetBField()*10;
385 for (int ib=0;ib<nbg;ib++) {
386 trBg = *((AliExternalTrackParam*)bgPoolITS->At(ib));
387 trBg.GetXYZ(posBg); // lab position
390 alpBg = trBg.GetAlpha();
391 parBg = (double*)trBg.GetParameter();
392 covBg = (double*)trBg.GetCovariance();
394 memcpy(&parBgOr,parBg,5*sizeof(double));
395 memcpy(&covBgOr,covBg,15*sizeof(double));
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);
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);
411 trBg.Set(xBg,alpBg,parBgOr, covBgOr); // restore original track
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);
423 hchi->Scale(nCorr/nTest);
428 //_________________________________________________________________________________________
429 TH1F* PrepChi2Ndf(int ndf, int ngen)
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);
439 //____________________________________________________
440 Double_t CalcCorrProb(TH1* chiSig, TH1* chiBg, double &smBgRet, double maxChi)
442 // Calculate probability of all fake matches being worse than correct one
443 // The chi2 histos should be normalized
444 int nb = chiSig->GetNbinsX();
446 int nt = chiSig->FindBin(maxChi);
450 double bgchiC = maxChi;
451 if (bgchiC<0) bgchiC = 5.;
452 int maxBgBinCut = chiSig->FindBin(bgchiC);
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;
463 smBg = 1.-TMath::Exp(smBg); // prob. to find bg within a cut in absence of correct match
467 //__________________________________________________
468 Bool_t PassFastCheck(AliExternalTrackParam &tr1, AliExternalTrackParam &tr2, double cut)
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};
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();
481 for (int i=0;i<5;i++) {
482 double df = par1[i]-par2[i];
483 double err2 = cov1[erId[i]]+cov2[erId[i]];
486 return (chi2>cut*kNDF) ? kFALSE : kTRUE;