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