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