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