1 #include "KMCDetector.h"
\r
4 #include <TProfile.h>
\r
6 #include <TMatrixD.h>
\r
9 #include <TFormula.h>
\r
10 #include <TCanvas.h>
\r
11 #include <TEllipse.h>
\r
13 #include <TRandom3.h>
\r
14 #include <TGraphErrors.h>
\r
15 #include <TStopwatch.h>
\r
17 /***********************************************************
\r
19 Fast Simulation tool for Inner Tracker Systems
\r
21 ***********************************************************/
\r
24 #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
\r
26 #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
\r
27 #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm)
\r
28 #define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)
\r
29 // #define dNdEtaCent 2300//15000 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known)
\r
31 #define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
\r
32 #define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
\r
33 #define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
\r
34 #define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0)
\r
35 #define EfficiencySearchFlag 2 // Define search method:
\r
36 // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
\r
38 #define PionMass 0.139 // Mass of the Pion
\r
39 #define KaonMass 0.498 // Mass of the Kaon
\r
40 #define D0Mass 1.865 // Mass of the D0
\r
42 //TMatrixD *probKomb; // table for efficiency kombinatorics
\r
46 Int_t KMCProbe::fgNITSLayers = 0;
\r
47 Double_t KMCProbe::fgMissingHitPenalty = 2.;
\r
48 //__________________________________________________________________________
\r
49 KMCProbe& KMCProbe::operator=(const KMCProbe& src)
\r
52 AliExternalTrackParam::operator=(src);
\r
56 fFakes = src.fFakes;
\r
57 fNHits = src.fNHits;
\r
58 fNHitsITS = src.fNHitsITS;
\r
59 fNHitsITSFake = src.fNHitsITSFake;
\r
60 fInnLrCheck = src.fInnLrCheck;
\r
61 for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
\r
66 //__________________________________________________________________________
\r
67 KMCProbe::KMCProbe(KMCProbe& src)
\r
68 : AliExternalTrackParam(src),
\r
74 fNHitsITS(src.fNHitsITS),
\r
75 fNHitsITSFake(src.fNHitsITSFake),
\r
76 fInnLrCheck(src.fInnLrCheck)
\r
78 for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
\r
81 //__________________________________________________________________________
\r
82 void KMCProbe::ResetCovMat()
\r
85 double *trCov = (double*)GetCovariance();
\r
86 double *trPars = (double*)GetParameter();
\r
87 const double kLargeErr2Coord = 50*50;
\r
88 const double kLargeErr2Dir = 0.6*0.6;
\r
89 const double kLargeErr2PtI = 0.5*0.5;
\r
90 for (int ic=15;ic--;) trCov[ic] = 0.;
\r
91 trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
\r
92 trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
\r
93 trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
\r
97 //__________________________________________________________________________
\r
98 void KMCProbe::Print(Option_t* option) const
\r
100 printf("M=%.3f Chi2=%7.2f (Norm:%6.2f|%d) Hits: Total:%d ITS:%d ITSFakes:%d | Y:%+8.4f Z: %+8.4f |",
\r
101 fMass,fChi2,GetNormChi2(kTRUE),fInnLrCheck,fNHits,fNHitsITS,fNHitsITSFake, GetY(),GetZ());
\r
102 for (int i=0;i<fgNITSLayers;i++) {
\r
103 if (!IsHit(i)) printf(".");
\r
104 else printf("%c",IsHitFake(i) ? '-':'+');
\r
106 printf("|%s\n",IsKilled() ? " KILLED":"");
\r
107 TString opt = option;
\r
108 if (!opt.IsNull()) AliExternalTrackParam::Print(option);
\r
111 //__________________________________________________________________________
\r
112 Int_t KMCProbe::Compare(const TObject* obj) const
\r
114 // compare to tracks
\r
115 const KMCProbe* trc = (KMCProbe*) obj;
\r
116 if (trc->IsKilled()) {
\r
117 if (IsKilled()) return 0;
\r
120 else if (IsKilled()) return 1;
\r
121 double chi2a = GetNormChi2(kTRUE);
\r
122 double chi2b = trc->GetNormChi2(kTRUE);
\r
123 if (chi2a<chi2b) return -1;
\r
124 if (chi2a>chi2b) return 1;
\r
128 //__________________________________________________________________________
\r
129 Bool_t KMCProbe::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
\r
131 // Get local X of the track position estimated at the radius lab radius r.
\r
132 // The track curvature is accounted exactly
\r
134 // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
\r
135 // 0 - take the intersection closest to the current track position
\r
136 // >0 - go along the track (increasing fX)
\r
137 // <0 - go backward (decreasing fX)
\r
139 // special case of R=0
\r
140 if (r<kAlmost0) {x=0; return kTRUE;}
\r
142 const double* pars = GetParameter();
\r
143 const Double_t &fy=pars[0], &sn = pars[2];
\r
145 double fx = GetX();
\r
146 double crv = GetC(bz);
\r
147 if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
\r
148 if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
\r
149 double det = (r-fx)*(r+fx);
\r
150 if (det<0) return kFALSE; // does not reach raduis r
\r
152 if (dir==0) return kTRUE;
\r
153 det = TMath::Sqrt(det);
\r
154 if (dir>0) { // along the track direction
\r
155 if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
\r
156 else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
\r
158 else if(dir>0) { // agains track direction
\r
159 if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
\r
160 else if (fy>det) return kFALSE; // track is against Y axis
\r
163 else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
\r
164 double det = (r-fy)*(r+fy);
\r
165 if (det<0) return kFALSE; // does not reach raduis r
\r
166 det = TMath::Sqrt(det);
\r
168 x = fx>0 ? det : -det; // choose the solution requiring the smalest step
\r
171 else if (dir>0) { // along the track direction
\r
172 if (fx > det) return kFALSE; // current point is in on the right from the circle
\r
173 else if (fx <-det) x = -det; // on the left
\r
174 else x = det; // within the circle
\r
176 else { // against the track direction
\r
177 if (fx <-det) return kFALSE;
\r
178 else if (fx > det) x = det;
\r
182 else { // general case of straight line
\r
183 double cs = TMath::Sqrt((1-sn)*(1+sn));
\r
184 double xsyc = fx*sn-fy*cs;
\r
185 double det = (r-xsyc)*(r+xsyc);
\r
186 if (det<0) return kFALSE; // does not reach raduis r
\r
187 det = TMath::Sqrt(det);
\r
188 double xcys = fx*cs+fy*sn;
\r
190 if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
\r
191 else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
\r
192 if (t>=-det) t += -det; // take minimal step giving t>0
\r
193 else return kFALSE; // both solutions have negative t
\r
195 else { // go in increasing fx direction. (t+-det < 0)
\r
196 if (t<det) t -= det; // take minimal step giving t<0
\r
197 else return kFALSE; // both solutions have positive t
\r
203 // get center of the track circle
\r
204 double tR = 1./crv; // track radius (for the moment signed)
\r
205 double cs = TMath::Sqrt((1-sn)*(1+sn));
\r
206 double x0 = fx - sn*tR;
\r
207 double y0 = fy + cs*tR;
\r
208 double r0 = TMath::Sqrt(x0*x0+y0*y0);
\r
209 // printf("Xc:%+e Yc:%+e\n",x0,y0);
\r
211 if (r0<=kAlmost0) {
\r
212 AliDebug(2,Form("r0 = %f",r0));
\r
214 } // the track is concentric to circle
\r
215 tR = TMath::Abs(tR);
\r
216 double tR2r0 = tR/r0;
\r
217 double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
\r
218 double det = (1.-g)*(1.+g);
\r
220 AliDebug(2,Form("g=%f tR=%f r0=%f\n",g,tR, r0));
\r
222 } // does not reach raduis r
\r
223 det = TMath::Sqrt(det);
\r
225 // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
\r
226 // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
\r
227 // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
\r
229 double tmp = 1.+g*tR2r0;
\r
232 if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
\r
233 double dfx = tR2r0*TMath::Abs(y0)*det;
\r
234 double dfy = tR2r0*x0*TMath::Sign(det,y0);
\r
235 if (dir==0) { // chose the one which corresponds to smallest step
\r
236 double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
\r
237 if (delta<0) x += dfx;
\r
240 else if (dir>0) { // along track direction: x must be > fx
\r
241 x -= dfx; // try the smallest step (dfx is positive)
\r
242 if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;
\r
244 else { // backward: x must be < fx
\r
245 x += dfx; // try the smallest step (dfx is positive)
\r
246 if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;
\r
249 else { // special case: track touching the circle just in 1 point
\r
250 if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE;
\r
257 //____________________________________
\r
258 Bool_t KMCProbe::PropagateToR(double r, double b, int dir)
\r
265 const double kTiny = 1e-4;
\r
267 if (!GetXatLabR(r ,xR, b, dir)) {
\r
268 // printf("Track with pt=%f cannot reach radius %f\n",Pt(),r);
\r
273 if (!PropagateTo(xR, b)) {
\r
274 if (AliLog::GetGlobalDebugLevel()>2) {
\r
275 printf("Failed to propagate to X=%f for R=%f\n",xR,r);
\r
280 double rcurr2 = xR*xR + GetY()*GetY();
\r
281 if (TMath::Abs(rcurr2-rr)<kTiny || rr<kAlmost0) return kTRUE;
\r
283 // two radii correspond to this X...
\r
284 double pos[3]; GetXYZ(pos);
\r
285 double phi = TMath::ATan2(pos[1],pos[0]);
\r
286 if (!Rotate(phi)) {
\r
287 if (AliLog::GetGlobalDebugLevel()>2) {
\r
288 printf("Failed to rotate to %f to propagate to R=%f\n",phi,r);
\r
294 if (AliLog::GetGlobalDebugLevel()>2) {
\r
295 printf("Failed to propagate to R=%f after %d steps\n",r,iter);
\r
305 //__________________________________________________________________________
\r
306 Bool_t KMCProbe::CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward)
\r
308 // printf("before at r=%.1f p=%.4f\n",lr->fR, P());
\r
309 if (AliExternalTrackParam::CorrectForMeanMaterial(lr->fx2X0, inward ? lr->fXRho : -lr->fXRho, GetMass() , kTRUE)) {
\r
310 // printf("after at r=%.1f p=%.4f\n",lr->fR, P());
\r
313 AliDebug(2,Form("Failed to apply material correction, X/X0=%.4f", lr->fx2X0));
\r
314 if (AliLog::GetGlobalDebugLevel()>1) Print();
\r
318 /////////////////////////////////////////////////////////////////////////////
\r
319 ClassImp(KMCCluster)
\r
321 //_________________________________________________________________________
\r
322 KMCCluster::KMCCluster(KMCCluster &src)
\r
324 fY(src.fY),fZ(src.fZ),fX(src.fX),fPhi(src.fPhi)
\r
327 //__________________________________________________________________________
\r
328 KMCCluster& KMCCluster::operator=(const KMCCluster& src)
\r
331 TObject::operator=(src);
\r
340 //_________________________________________________________________________
\r
341 void KMCCluster::Print(Option_t *) const
\r
343 printf(" Local YZ = (%3.4lf,%3.4lf) | X=%3.4lf phi: %+.3f %s\n",fY,fZ,fX,fPhi,IsKilled()?"Killed":"");
\r
346 /////////////////////////////////////////////////////////////////////////////
\r
349 Double_t KMCLayer::fgDefEff = 1.0;
\r
350 //__________________________________________________________________________
\r
351 KMCLayer::KMCLayer(char *name) :
\r
352 TNamed(name,name),fR(0),fx2X0(0),fPhiRes(0),fZRes(0),fEff(0),fIsDead(kFALSE),fType(-1),fActiveID(-1),fSig2EstD(999),fSig2EstZ(999),
\r
353 fClCorr(),fClMC(),fClBg("KMCCluster",5), fTrCorr(), fTrMC("KMCProbe",5)
\r
358 //__________________________________________________________________________
\r
359 void KMCLayer::Reset()
\r
364 fSig2EstD = fSig2EstZ = 999;
\r
368 //__________________________________________________________________________
\r
369 KMCProbe* KMCLayer::AddMCTrack(KMCProbe* src)
\r
371 int ntr = GetNMCTracks();
\r
373 if (src) prb = new(fTrMC[ntr]) KMCProbe(*src);
\r
374 else prb = new(fTrMC[ntr]) KMCProbe();
\r
375 if (!IsDead()) prb->ResetHit(GetActiveID());
\r
379 //__________________________________________________________________________
\r
380 void KMCLayer::Print(Option_t *opt) const
\r
382 printf("Lr%3d(A%3d) %10s R=%5.1f X2X0=%.3f XRho=%.3f SigY=%.4f SigZ=%.4f Eff:%4.2f\n",GetUniqueID(),fActiveID,GetName(), fR, fx2X0,fXRho,fPhiRes,fZRes,fEff);
\r
383 TString opts = opt; opts.ToLower();
\r
384 if (opts.Contains("c")) {
\r
385 printf("Clusters: MC: %+7.4f:%+7.4f Ideal: %+7.4f:%+7.4f NBgCl: %3d NTrMC: %4d\n",fClMC.fY,fClMC.fZ, fClCorr.fY,fClCorr.fZ, GetNBgClusters(),GetNMCTracks());
\r
389 /////////////////////////////////////////////////////////////////////////////
\r
390 Double_t KMCDetector::fgVtxConstraint[2]={-1,-1};
\r
392 ClassImp(KMCDetector)
\r
393 KMCDetector::KMCDetector() :
\r
394 TNamed("test_detector","detector"),
\r
397 fNActiveITSLayers(0),
\r
398 fLastActiveLayer(-1),
\r
399 fLastActiveITSLayer(-1),
\r
400 fLastActiveLayerTracked(-1),
\r
402 fLhcUPCscale(1.0),
\r
403 fIntegrationTime(0.02), // in ms
\r
404 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
\r
405 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
\r
406 fParticleMass(0.140), // Standard: pion mass
\r
407 fMaxRadiusSlowDet(10.),
\r
408 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
\r
409 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
\r
412 fMaxNormChi2NDF(5.),
\r
415 fMaxChi2ClSQ(4.), // precalulated internally
\r
416 fMaxSeedToPropagate(300),
\r
418 fUseBackground(kFALSE),
\r
419 fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
\r
420 fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
\r
421 fDensFactorEta(1.),
\r
424 fHMCLrResidRPhi(0),
\r
430 // default constructor
\r
432 // fLayers = new TObjArray();
\r
433 RequireMaxChi2Cl(fMaxChi2Cl); // just to precalulate default square
\r
436 KMCDetector::KMCDetector(char *name, char *title)
\r
437 : TNamed(name,title),
\r
440 fNActiveITSLayers(0),
\r
441 fLastActiveLayer(-1),
\r
442 fLastActiveITSLayer(-1),
\r
443 fLastActiveLayerTracked(-1),
\r
446 fIntegrationTime(0.02), // in ms
\r
447 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
\r
448 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
\r
449 fParticleMass(0.140), // Standard: pion mass
\r
450 fMaxRadiusSlowDet(10.),
\r
451 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
\r
452 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
\r
455 fMaxNormChi2NDF(5.),
\r
458 fMaxChi2ClSQ(3.), // precalulated internally
\r
459 fMaxSeedToPropagate(50),
\r
461 fUseBackground(kFALSE),
\r
462 fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
\r
463 fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
\r
464 fDensFactorEta(1.),
\r
467 fHMCLrResidRPhi(0),
\r
473 // default constructor, that set the name and title
\r
475 // fLayers = new TObjArray();
\r
477 KMCDetector::~KMCDetector() { //
\r
478 // virtual destructor
\r
483 void KMCDetector::InitMCWatchHistos()
\r
485 // init utility histos used for MC tuning
\r
486 enum {kNBinRes=1000};
\r
487 const double kMaxResidRPhi=1.0,kMaxResidZ=1.0,kMaxChi2=50;
\r
488 int nlr = fNActiveITSLayers;
\r
489 TString nam = "mc_residrhi";
\r
490 TString tit = "MC $Delta Cl-Tr R#phi";
\r
491 fHMCLrResidRPhi = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidRPhi,nlr,-0.5,nlr-0.5);
\r
492 fHMCLrResidRPhi->GetXaxis()->SetTitle("cl-tr #Delta r#phi");
\r
493 fHMCLrResidRPhi->Sumw2();
\r
496 tit = "MC $Delta Cl-Tr Z";
\r
497 fHMCLrResidZ = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidZ,nlr,-0.5,nlr-0.5);
\r
498 fHMCLrResidZ->GetXaxis()->SetTitle("cl-tr #Delta z");
\r
499 fHMCLrResidZ->Sumw2();
\r
502 tit = "MC #chi^{2} Cl-Tr Z";
\r
503 fHMCLrChi2 = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxChi2,nlr,-0.5,nlr-0.5);
\r
504 fHMCLrChi2->GetXaxis()->SetTitle("cl-tr #chi^{2}");
\r
505 fHMCLrChi2->Sumw2();
\r
507 SetBit(kUtilHisto);
\r
511 void KMCDetector::AddLayer(char *name, Float_t radius, Float_t x2X0, Float_t xrho, Float_t phiRes, Float_t zRes, Float_t eff) {
\r
513 // Add additional layer to the list of layers (ordered by radius)
\r
516 KMCLayer *newLayer = (KMCLayer*) fLayers.FindObject(name);
\r
519 newLayer = new KMCLayer(name);
\r
520 newLayer->fR = radius;
\r
521 newLayer->fx2X0 = x2X0;
\r
522 newLayer->fXRho = xrho;
\r
523 newLayer->fPhiRes = phiRes;
\r
524 newLayer->fZRes = zRes;
\r
525 eff = TMath::Min(1.f,eff);
\r
526 newLayer->fEff = eff <0 ? KMCLayer::GetDefEff() : eff;
\r
527 newLayer->fActiveID = -2;
\r
528 TString lname = name;
\r
529 newLayer->fType = KMCLayer::kTypeNA;
\r
530 if (lname.Contains("tpc")) newLayer->fType = KMCLayer::kTPC;
\r
531 else if (lname.Contains("its")) newLayer->fType = KMCLayer::kITS;
\r
532 if (lname.Contains("vertex")) newLayer->SetBit(KMCLayer::kBitVertex);
\r
534 if (newLayer->fType==KMCLayer::kTypeNA) printf("Attention: the layer %s has undefined type\n",name);
\r
536 newLayer->fIsDead = (newLayer->fPhiRes==RIDICULOUS && newLayer->fZRes==RIDICULOUS);
\r
538 if (fLayers.GetEntries()==0)
\r
539 fLayers.Add(newLayer);
\r
542 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
543 KMCLayer *l = (KMCLayer*)fLayers.At(i);
\r
544 if (radius<l->fR) { fLayers.AddBefore(l,newLayer); break; }
\r
545 if (radius>l->fR && (i+1)==fLayers.GetEntries() ) fLayers.Add(newLayer); // even bigger then last one
\r
553 printf("Layer with the name %s does already exist\n",name);
\r
557 //____________________________________________________________
\r
558 void KMCDetector::ClassifyLayers()
\r
560 // assign active Id's, etc
\r
561 fLastActiveLayer = -1;
\r
562 fLastActiveITSLayer = -1;
\r
563 fNActiveLayers = 0;
\r
564 fNActiveITSLayers = 0;
\r
566 int nl = GetNLayers();
\r
567 for (int il=0;il<nl;il++) {
\r
568 KMCLayer* lr = GetLayer(il);
\r
569 lr->SetUniqueID(il);
\r
570 if (!lr->IsDead()) {
\r
571 fLastActiveLayer = il;
\r
572 lr->fActiveID = fNActiveLayers++;
\r
574 fLastActiveITSLayer = il;
\r
575 fNActiveITSLayers++;
\r
580 KMCProbe::SetNITSLayers(fNActiveITSLayers);
\r
583 //____________________________________________________________
\r
584 void KMCDetector::KillLayer(char *name) {
\r
586 // Marks layer as dead. Contribution only by Material Budget
\r
589 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
591 printf("Layer %s not found - cannot mark as dead\n",name);
\r
593 tmp->fPhiRes = 999999;
\r
594 tmp->fZRes = 999999;
\r
595 tmp->fIsDead = kTRUE;
\r
600 void KMCDetector::SetRadius(char *name, Float_t radius) {
\r
602 // Set layer radius [cm]
\r
604 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
607 printf("Layer %s not found - cannot set radius\n",name);
\r
610 Float_t tmpRadL = tmp->fx2X0;
\r
611 Float_t tmpPhiRes = tmp->fPhiRes;
\r
612 Float_t tmpZRes = tmp->fZRes;
\r
613 Float_t tmpXRho = tmp->fXRho;
\r
614 RemoveLayer(name); // so that the ordering is correct
\r
615 AddLayer(name,radius,tmpRadL,tmpXRho,tmpPhiRes,tmpZRes);
\r
619 Float_t KMCDetector::GetRadius(char *name) {
\r
621 // Return layer radius [cm]
\r
624 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
626 printf("Layer %s not found - cannot get radius\n",name);
\r
633 void KMCDetector::SetRadiationLength(char *name, Float_t x2X0) {
\r
635 // Set layer material [cm]
\r
638 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
640 printf("Layer %s not found - cannot set layer material\n",name);
\r
646 Float_t KMCDetector::GetRadiationLength(char *name) {
\r
648 // Return layer radius [cm]
\r
651 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
653 printf("Layer %s not found - cannot get layer material\n",name);
\r
661 void KMCDetector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
\r
663 // Set layer resolution in [cm]
\r
666 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
668 printf("Layer %s not found - cannot set resolution\n",name);
\r
670 tmp->fPhiRes = phiRes;
\r
672 tmp->fIsDead = (zRes==RIDICULOUS && phiRes==RIDICULOUS);
\r
677 Float_t KMCDetector::GetResolution(char *name, Int_t axis) {
\r
679 // Return layer resolution in [cm]
\r
680 // axis = 0: resolution in rphi
\r
681 // axis = 1: resolution in z
\r
684 KMCLayer *tmp = GetLayer(name);
\r
686 printf("Layer %s not found - cannot get resolution\n",name);
\r
688 if (axis==0) return tmp->fPhiRes;
\r
689 if (axis==1) return tmp->fZRes;
\r
690 printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
\r
695 void KMCDetector::SetLayerEfficiency(char *name, Float_t eff) {
\r
697 // Set layer efficnecy (prop that his is missed within this layer)
\r
700 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
702 printf("Layer %s not found - cannot set layer efficiency\n",name);
\r
708 Float_t KMCDetector::GetLayerEfficiency(char *name) {
\r
710 // Get layer efficnecy (prop that his is missed within this layer)
\r
713 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
715 printf("Layer %s not found - cannot get layer efficneicy\n",name);
\r
723 void KMCDetector::RemoveLayer(char *name) {
\r
725 // Removes a layer from the list
\r
728 KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
\r
730 printf("Layer %s not found - cannot remove it\n",name);
\r
732 fLayers.Remove(tmp);
\r
738 void KMCDetector::PrintLayout() {
\r
740 // Prints the detector layout
\r
743 printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
\r
745 if (fLayers.GetEntries()>0)
\r
746 printf(" Name \t\t r [cm] \t X0 \t phi & z res [um]\n");
\r
749 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
750 tmp = (KMCLayer*)fLayers.At(i);
\r
752 // don't print all the tpc layers
\r
753 TString name(tmp->GetName());
\r
754 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
\r
756 printf("%d. %s \t %03.2f \t%1.4f\t ",i,
\r
757 tmp->GetName(), tmp->fR, tmp->fx2X0);
\r
758 if (tmp->fPhiRes==RIDICULOUS)
\r
761 printf("%3.0f ",tmp->fPhiRes*10000);
\r
762 if (tmp->fZRes==RIDICULOUS)
\r
765 printf("%3.0f\n",tmp->fZRes*10000);
\r
769 void KMCDetector::PlotLayout(Int_t plotDead) {
\r
771 // Plots the detector layout in Front view
\r
774 Double_t x0=0, y0=0;
\r
776 TGraphErrors *gr = new TGraphErrors();
\r
777 gr->SetPoint(0,0,0);
\r
778 KMCLayer *lastLayer = (KMCLayer*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->fR;
\r
779 gr->SetPointError(0,maxRad,maxRad);
\r
784 for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
\r
785 tmp = (KMCLayer*)fLayers.At(i);
\r
788 Double_t txtpos = tmp->fR;
\r
789 if ((tmp->IsDead())) txtpos*=-1; //
\r
790 TText *txt = new TText(x0,txtpos,tmp->GetName());
\r
791 txt->SetTextSizePixels(5); txt->SetTextAlign(21);
\r
792 if (!tmp->IsDead() || plotDead) txt->Draw();
\r
794 TEllipse *layEl = new TEllipse(x0,y0,tmp->fR);
\r
795 // layEl->SetFillColor(5);
\r
796 layEl->SetFillStyle(5001);
\r
797 layEl->SetLineStyle(tmp->IsDead()+1); // dashed if not active
\r
798 layEl->SetLineColor(4);
\r
799 TString name(tmp->GetName());
\r
800 if (!tmp->IsDead()) layEl->SetLineWidth(2);
\r
801 if (name.Contains("tpc") ) layEl->SetLineColor(29);
\r
803 if (!tmp->IsDead() || plotDead) layEl->Draw();
\r
811 void KMCDetector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
\r
813 // Emulates the TPC
\r
815 // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
\r
818 AddLayer((char*)"IFCtpc", 77.8,0.01367, 0); // Inner Field cage (RS: set correct xrho for eloss)
\r
820 // % Radiation Lengths ... Average per TPC row (i.e. total/159 )
\r
821 Float_t x2X0PerRow = 0.000036;
\r
823 Float_t tpcInnerRadialPitch = 0.75 ; // cm
\r
824 Float_t tpcMiddleRadialPitch = 1.0 ; // cm
\r
825 Float_t tpcOuterRadialPitch = 1.5 ; // cm
\r
826 // Float_t tpcInnerPadWidth = 0.4 ; // cm
\r
827 // Float_t tpcMiddlePadWidth = 0.6 ; // cm
\r
828 // Float_t tpcOuterPadWidth = 0.6 ; // cm
\r
829 Float_t innerRows = 63 ;
\r
830 Float_t middleRows = 64 ;
\r
831 Float_t outerRows = 32 ;
\r
832 Float_t tpcRows = (innerRows + middleRows + outerRows) ;
\r
833 Float_t rowOneRadius = 85.2 ; // cm
\r
834 Float_t row64Radius = 135.1 ; // cm
\r
835 Float_t row128Radius = 199.2 ; // cm
\r
837 for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
\r
839 Float_t rowRadius =0;
\r
841 rowRadius = rowOneRadius + k*tpcInnerRadialPitch ;
\r
842 else if ( k>=innerRows && k<(innerRows+middleRows) )
\r
843 rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
\r
844 else if (k>=(innerRows+middleRows) && k<tpcRows )
\r
845 rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
\r
848 AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0, phiResMean,zResMean);
\r
850 AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0); // non "active" row
\r
857 void KMCDetector::RemoveTPC() {
\r
859 // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
\r
861 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
862 tmp = (KMCLayer*)fLayers.At(i);
\r
863 TString name(tmp->GetName());
\r
864 if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
\r
866 RemoveLayer((char*)"IFC");
\r
871 Double_t KMCDetector::ThetaMCS ( Double_t mass, Double_t x2X0, Double_t momentum ) const
\r
874 // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
\r
877 Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ;
\r
878 Double_t theta = 0.0 ; // Momentum and mass in GeV
\r
879 // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
\r
880 if ( x2X0 > 0 ) theta = 0.0136 * TMath::Sqrt(x2X0) / ( beta * momentum ) * (1+0.038*TMath::Log(x2X0)) ;
\r
885 Double_t KMCDetector::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
887 // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
\r
888 // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
\r
889 // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
\r
890 Double_t sx, sy, goodHit ;
\r
891 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
\r
892 sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ;
\r
893 goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ;
\r
894 return ( goodHit ) ;
\r
898 Double_t KMCDetector::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
900 // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
\r
901 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
\r
902 Double_t sx, goodHit ;
\r
903 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
\r
904 goodHit = 1./(1+sx) ;
\r
905 return ( goodHit ) ;
\r
908 Double_t KMCDetector::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
910 // Based on work by Ruben Shahoyen
\r
911 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
\r
912 // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
\r
913 // Following is correct for 2 DOF
\r
915 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
\r
916 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
\r
917 Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));
\r
918 return ( goodHit ) ;
\r
921 Double_t KMCDetector::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
923 // Based on work by Ruben Shahoyen
\r
924 // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)
\r
926 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
\r
927 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
\r
928 Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));
\r
929 return ( nullHit ) ;
\r
932 Double_t KMCDetector::HitDensity ( Double_t radius )
\r
934 // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
\r
935 // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
\r
936 // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
\r
937 // Note that this function assumes we are working in CM and CM**2 [not meters].
\r
938 // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
\r
940 // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors
\r
941 // and only fast detectors reside outside this radius.
\r
942 Double_t arealDensity = 0 ;
\r
943 if (radius<0.01) return 0;
\r
945 if ( radius >= fMaxRadiusSlowDet )
\r
947 arealDensity = OneEventHitDensity(fdNdEtaCent,radius) ; // Fast detectors see central collision density (only)
\r
948 arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background
\r
951 if (radius < fMaxRadiusSlowDet )
\r
952 { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
\r
953 arealDensity = OneEventHitDensity(fdNdEtaCent,radius)
\r
954 + IntegratedHitDensity(dNdEtaMinB,radius)
\r
955 + UpcHitDensity(radius) ;
\r
956 arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
\r
957 // Increase density due to background
\r
960 return ( arealDensity ) ;
\r
964 double KMCDetector::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
\r
966 // This is for one event at the vertex. No smearing.
\r
967 double den = multiplicity / (2.*TMath::Pi()*radius*radius) * fDensFactorEta ; // 2 eta ?
\r
968 // note: surface of sphere is '4*pi*r^2'
\r
969 // surface of cylinder is '2*pi*r* h'
\r
974 double KMCDetector::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
\r
976 // The integral of minBias events smeared over a gaussian vertex distribution.
\r
977 // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
\r
979 Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
\r
980 Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
\r
982 // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
\r
983 double dens1 = OneEventHitDensity(multiplicity,radius);
\r
984 if ( den < dens1 ) den = dens1;
\r
990 double KMCDetector::UpcHitDensity(Double_t radius)
\r
992 // QED electrons ...
\r
994 Double_t mUPCelectrons ; ;
\r
995 // mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
\r
996 mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB; // Fit to 'Rossegger,Sadovsky'-Alice simulation
\r
997 if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R
\r
998 mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate
\r
999 mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC
\r
1001 return mUPCelectrons ;
\r
1005 double KMCDetector::Dist(double z, double r)
\r
1007 // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution
\r
1008 // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
\r
1009 // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
\r
1010 Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration
\r
1011 Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work
\r
1012 double dist = 0.0 ;
\r
1013 double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm
\r
1014 double z0 = 0.0 ; //cm
\r
1015 for(int i=0; i<nsteps; i++){
\r
1016 if ( i == nsteps-1 ) index = 1 ;
\r
1017 z0 = -4*SigmaD + i*dz0 ;
\r
1018 dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
\r
1019 (1/sqrt((z-z0)*(z-z0) + r*r)) ;
\r
1020 if ( index != 4 ) index = 4; else index = 2 ;
\r
1025 #define PZero 0.861 // Momentum of back to back decay particles in the CM frame
\r
1026 #define EPiZero 0.872 // Energy of the pion from a D0 decay at rest
\r
1027 #define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest
\r
1029 Double_t KMCDetector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) const {
\r
1030 // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi.
\r
1032 Double_t const1 = pt / D0Mass ;
\r
1033 Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
\r
1034 Double_t sum, ptPi, ptK ;
\r
1035 Double_t effp, effk ;
\r
1038 for ( Int_t k = 0 ; k < 360 ; k++ ) {
\r
1040 Double_t theta = k * TMath::Pi() / 180. ;
\r
1042 ptPi = TMath::Sqrt(
\r
1043 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
\r
1044 const1*const1*EPiZero*EPiZero -
\r
1045 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
\r
1046 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
\r
1049 ptK = TMath::Sqrt(
\r
1050 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
\r
1051 const1*const1*EKZero*EKZero +
\r
1052 2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
\r
1053 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
\r
1056 // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
\r
1057 Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ;
\r
1058 Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ;
\r
1060 if ( pionindex >= 20 ) pionindex = 399 ;
\r
1061 if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ;
\r
1062 if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
\r
1063 if ( effp < 0 ) effp = 0 ;
\r
1065 if ( kaonindex >= 20 ) kaonindex = 399 ;
\r
1066 if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ;
\r
1067 if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
\r
1068 if ( effk < 0 ) effk = 0 ;
\r
1070 // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
\r
1072 sum += effp * effk ;
\r
1076 Double_t mean =sum/360;
\r
1081 KMCProbe* KMCDetector::PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi, double x,double y, double z)
\r
1083 // Prepare trackable Kalman track at the farthest position
\r
1085 // Set track parameters
\r
1086 // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
\r
1088 fProbe.SetMass(mass);
\r
1089 KMCProbe* probe = new KMCProbe(fProbe);
\r
1090 double *trPars = (double*)probe->GetParameter();
\r
1091 double *trCov = (double*)probe->GetCovariance();
\r
1092 double xyz[3] = {x,y,z};
\r
1093 probe->Global2LocalPosition(xyz,phi);
\r
1094 probe->Set(xyz[0],phi,trPars,trCov);
\r
1095 trPars[KMCProbe::kY] = xyz[1];
\r
1096 trPars[KMCProbe::kZ] = xyz[2];
\r
1097 trPars[KMCProbe::kSnp] = 0; // track along X axis at the vertex
\r
1098 trPars[KMCProbe::kTgl] = TMath::Tan(lambda); // dip
\r
1099 trPars[KMCProbe::kPtI] = charge/pt; // q/pt
\r
1101 // put tiny errors to propagate to the outer-most radius
\r
1102 trCov[KMCProbe::kY2] = trCov[KMCProbe::kZ2] = trCov[KMCProbe::kSnp2] = trCov[KMCProbe::kTgl2] = trCov[KMCProbe::kPtI2] = 1e-20;
\r
1103 fProbe = *probe; // store original track
\r
1105 // propagate to last layer
\r
1106 fLastActiveLayerTracked = 0;
\r
1107 for (Int_t j=0; j<=fLastActiveLayer; j++) {
\r
1108 KMCLayer* lr = GetLayer(j);
\r
1111 if (!PropagateToLayer(probe,lr,1)) break;
\r
1112 if (!probe->CorrectForMeanMaterial(lr, kFALSE)) break;
\r
1114 lr->fClCorr.Set(probe->GetY(),probe->GetZ(), probe->GetX(), probe->GetAlpha());
\r
1115 if (!lr->IsDead()) fLastActiveLayerTracked = j;
\r
1117 probe->ResetCovMat();// reset cov.matrix
\r
1118 printf("Last active layer trracked: %d (out of %d)\n",fLastActiveLayerTracked,fLastActiveLayer);
\r
1125 TGraph * KMCDetector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
\r
1127 // returns the momentum resolution
\r
1130 TGraph *graph = new TGraph(20, fTransMomenta, fMomentumRes);
\r
1131 graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
\r
1132 // graph->GetXaxis()->SetRangeUser(0.,5.0) ;
\r
1133 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1134 graph->GetXaxis()->CenterTitle();
\r
1135 graph->GetXaxis()->SetNoExponent(1) ;
\r
1136 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1137 graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;
\r
1138 graph->GetYaxis()->CenterTitle();
\r
1140 graph->SetMaximum(20) ;
\r
1141 graph->SetMinimum(0.1) ;
\r
1142 graph->SetLineColor(color);
\r
1143 graph->SetMarkerColor(color);
\r
1144 graph->SetLineWidth(linewidth);
\r
1150 TGraph * KMCDetector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
\r
1152 // Returns the pointing resolution
\r
1153 // axis = 0 ... rphi pointing resolution
\r
1154 // axis = 1 ... z pointing resolution
\r
1157 TGraph * graph = 0;
\r
1160 graph = new TGraph ( 20, fTransMomenta, fResolutionRPhi ) ;
\r
1161 graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
\r
1162 graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
\r
1164 graph = new TGraph ( 20, fTransMomenta, fResolutionZ ) ;
\r
1165 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
\r
1166 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
\r
1169 graph->SetMinimum(1) ;
\r
1170 graph->SetMaximum(300.1) ;
\r
1171 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1172 graph->GetXaxis()->CenterTitle();
\r
1173 graph->GetXaxis()->SetNoExponent(1) ;
\r
1174 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1175 graph->GetYaxis()->CenterTitle();
\r
1177 graph->SetLineWidth(linewidth);
\r
1178 graph->SetLineColor(color);
\r
1179 graph->SetMarkerColor(color);
\r
1186 TGraph * KMCDetector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
\r
1188 // returns the Pointing resolution (accoring to Telescope equation)
\r
1189 // axis =0 ... in rphi
\r
1190 // axis =1 ... in z
\r
1193 Double_t resolution[20];
\r
1195 Double_t layerResolution[2];
\r
1196 Double_t layerRadius[2];
\r
1197 Double_t layerThickness[2];
\r
1199 Int_t count =0; // search two first active layers
\r
1200 printf("Telescope equation for layers: ");
\r
1201 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
1202 KMCLayer *l = (KMCLayer*)fLayers.At(i);
\r
1203 if (!l->IsDead() && l->fR>0) {
\r
1204 layerRadius[count] = l->fR;
\r
1205 layerThickness[count] = l->fx2X0;
\r
1207 layerResolution[count] = l->fPhiRes;
\r
1209 layerResolution[count] = l->fZRes;
\r
1211 printf("%s, ",l->GetName());
\r
1214 if (count>=2) break;
\r
1218 Double_t pt, momentum, thickness,aMCS ;
\r
1219 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1221 for ( Int_t i = 0 ; i < 20 ; i++ ) {
\r
1222 // Reference data as if first two layers were acting all alone
\r
1223 pt = fTransMomenta[i] ;
\r
1224 momentum = pt / TMath::Cos(lambda) ; // Total momentum
\r
1225 resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
\r
1226 + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
\r
1227 resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
\r
1228 thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
\r
1229 aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
\r
1230 resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
\r
1231 resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns
\r
1236 TGraph* graph = new TGraph ( 20, fTransMomenta, resolution ) ;
\r
1239 graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
\r
1240 graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;
\r
1242 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
\r
1243 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;
\r
1245 graph->SetMinimum(1) ;
\r
1246 graph->SetMaximum(300.1) ;
\r
1247 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1248 graph->GetXaxis()->CenterTitle();
\r
1249 graph->GetXaxis()->SetNoExponent(1) ;
\r
1250 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1251 graph->GetYaxis()->CenterTitle();
\r
1253 graph->SetLineColor(color);
\r
1254 graph->SetMarkerColor(color);
\r
1255 graph->SetLineStyle(kDashed);
\r
1256 graph->SetLineWidth(linewidth);
\r
1262 TGraph * KMCDetector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
\r
1264 // particle = 0 ... choosen particle (setted particleMass)
\r
1265 // particle = 1 ... Pion
\r
1266 // particle = 2 ... Kaon
\r
1267 // particle = 3 ... D0
\r
1269 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1271 Double_t particleEfficiency[20]; // with chosen particle mass
\r
1272 Double_t kaonEfficiency[20], pionEfficiency[20], d0efficiency[20];
\r
1273 Double_t partEfficiency[2][20];
\r
1275 if (particle != 0) {
\r
1276 // resulting Pion and Kaon efficiency scaled with overall efficiency
\r
1277 Double_t doNotDecayFactor;
\r
1278 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
\r
1280 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1281 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
\r
1282 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
\r
1283 if ( massloop == 1 ) { // KAON
\r
1284 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
\r
1285 kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
\r
1287 doNotDecayFactor = 1.0 ;
\r
1288 pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
\r
1290 partEfficiency[0][j] = pionEfficiency[j];
\r
1291 partEfficiency[1][j] = kaonEfficiency[j];
\r
1295 // resulting estimate of the D0 efficiency
\r
1296 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1297 d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
\r
1300 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1301 particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
\r
1302 // NOTE: Decay factor (see kaon) should be included to be realiable
\r
1306 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1307 pionEfficiency[j] *= 100;
\r
1308 kaonEfficiency[j] *= 100;
\r
1309 d0efficiency[j] *= 100;
\r
1310 particleEfficiency[j] *= 100;
\r
1313 TGraph * graph = 0;
\r
1314 if (particle==0) {
\r
1315 graph = new TGraph ( 20, fTransMomenta, particleEfficiency ) ; // choosen mass
\r
1316 graph->SetLineWidth(1);
\r
1317 } else if (particle==1) {
\r
1318 graph = new TGraph ( 20, fTransMomenta, pionEfficiency ) ;
\r
1319 graph->SetLineWidth(1);
\r
1320 } else if (particle ==2) {
\r
1321 graph = new TGraph ( 20, fTransMomenta, kaonEfficiency ) ;
\r
1322 graph->SetLineWidth(1);
\r
1323 } else if (particle ==3) {
\r
1324 graph = new TGraph ( 20, fTransMomenta, d0efficiency ) ;
\r
1325 graph->SetLineStyle(kDashed);
\r
1329 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1330 graph->GetXaxis()->CenterTitle();
\r
1331 graph->GetXaxis()->SetNoExponent(1) ;
\r
1332 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1333 graph->GetYaxis()->SetTitle("Efficiency (%)") ;
\r
1334 graph->GetYaxis()->CenterTitle();
\r
1336 graph->SetMinimum(0.01) ;
\r
1337 graph->SetMaximum(100) ;
\r
1339 graph->SetLineColor(color);
\r
1340 graph->SetMarkerColor(color);
\r
1341 graph->SetLineWidth(linewidth);
\r
1346 TGraph * KMCDetector::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {
\r
1348 // particle = 0 ... choosen particle (setted particleMass)
\r
1349 // particle = 1 ... Pion
\r
1350 // particle = 2 ... Kaon
\r
1353 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1355 Double_t particleFake[20]; // with chosen particle mass
\r
1356 Double_t kaonFake[20], pionFake[20];
\r
1357 Double_t partFake[2][20];
\r
1359 if (particle != 0) {
\r
1360 // resulting Pion and Kaon efficiency scaled with overall efficiency
\r
1361 Double_t doNotDecayFactor;
\r
1362 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
\r
1364 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1365 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
\r
1366 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
\r
1367 if ( massloop == 1 ) { // KAON
\r
1368 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
\r
1369 kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;
\r
1371 pionFake[j] = fFake[0][j] ;
\r
1373 partFake[0][j] = pionFake[j];
\r
1374 partFake[1][j] = kaonFake[j];
\r
1379 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1380 particleFake[j] = fFake[2][j];
\r
1381 // NOTE: Decay factor (see kaon) should be included to be realiable
\r
1385 for ( Int_t j = 0 ; j < 20 ; j++ ) {
\r
1386 pionFake[j] *= 100;
\r
1387 kaonFake[j] *= 100;
\r
1388 particleFake[j] *= 100;
\r
1391 TGraph * graph = 0;
\r
1392 if (particle==0) {
\r
1393 graph = new TGraph ( 20, fTransMomenta, particleFake ) ; // choosen mass
\r
1394 graph->SetLineWidth(1);
\r
1395 } else if (particle==1) {
\r
1396 graph = new TGraph ( 20, fTransMomenta, pionFake ) ;
\r
1397 graph->SetLineWidth(1);
\r
1398 } else if (particle ==2) {
\r
1399 graph = new TGraph ( 20, fTransMomenta, kaonFake ) ;
\r
1400 graph->SetLineWidth(1);
\r
1403 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1404 graph->GetXaxis()->CenterTitle();
\r
1405 graph->GetXaxis()->SetNoExponent(1) ;
\r
1406 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1407 graph->GetYaxis()->SetTitle("Fake (%)") ;
\r
1408 graph->GetYaxis()->CenterTitle();
\r
1410 graph->SetMinimum(0.01) ;
\r
1411 graph->SetMaximum(100) ;
\r
1413 graph->SetLineColor(color);
\r
1414 graph->SetMarkerColor(color);
\r
1415 graph->SetLineWidth(linewidth);
\r
1421 TGraph* KMCDetector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
\r
1423 // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
\r
1424 // mode 0: impact parameter (convolution of pointing and vertex resolution)
\r
1425 // mode 1: pointing resolution
\r
1426 // mode 2: vtx resolution
\r
1429 TGraph *graph = new TGraph();
\r
1431 // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
\r
1432 TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
\r
1433 TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
\r
1435 TGraph *trackRes = GetGraphPointingResolution(axis,1);
\r
1436 Double_t *pt = trackRes->GetX();
\r
1437 Double_t *trRes = trackRes->GetY();
\r
1438 for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
\r
1439 Double_t vtxRes = 0;
\r
1441 vtxRes = vtxResRPhi.Eval(pt[ip]);
\r
1443 vtxRes = vtxResZ.Eval(pt[ip]);
\r
1446 graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
\r
1447 else if (mode ==1)
\r
1448 graph->SetPoint(ip,pt[ip],trRes[ip]);
\r
1450 graph->SetPoint(ip,pt[ip],vtxRes);
\r
1453 graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
\r
1454 graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;
\r
1456 graph->SetMinimum(1) ;
\r
1457 graph->SetMaximum(300.1) ;
\r
1458 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1459 graph->GetXaxis()->CenterTitle();
\r
1460 graph->GetXaxis()->SetNoExponent(1) ;
\r
1461 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1462 graph->GetYaxis()->CenterTitle();
\r
1464 graph->SetLineColor(color);
\r
1465 graph->SetMarkerColor(color);
\r
1466 graph->SetLineWidth(linewidth);
\r
1472 TGraph* KMCDetector::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
\r
1474 // returns graph according to the number
\r
1478 return GetGraphPointingResolution(0,color, linewidth); // dr
\r
1480 return GetGraphPointingResolution(1,color, linewidth); // dz
\r
1482 return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
\r
1484 return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
\r
1486 return GetGraphMomentumResolution(color, linewidth); // pt resolution
\r
1488 return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle
\r
1490 return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion
\r
1492 return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon
\r
1494 return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0
\r
1496 return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle
\r
1498 return GetGraphRecoFakes(1, color, linewidth); // Fake pion
\r
1500 return GetGraphRecoFakes(2, color, linewidth); // Fake kaon
\r
1502 printf(" Error: chosen graph number not valid\n");
\r
1508 void KMCDetector::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon, int setVer) {
\r
1510 // All New configuration with X0 = 0.3 and resolution = 4 microns
\r
1512 AddLayer((char*)"bpipe_its",2.0,0.0022, 0.092); // beam pipe, 0.5 mm Be
\r
1513 AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation
\r
1514 if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
\r
1515 printf("vertex will work as constraint: %.4f %.4f\n",fgVtxConstraint[0],fgVtxConstraint[1]);
\r
1516 SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
\r
1519 // new ideal Pixel properties?
\r
1520 Double_t x0 = 0.0050;
\r
1521 Double_t resRPhi = 0.0006;
\r
1522 Double_t resZ = 0.0006;
\r
1523 Double_t xrho = 0.0116; // assume 0.5mm of silicon for eloss
\r
1535 // default all new
\r
1537 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1538 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1539 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1540 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1541 AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
\r
1542 AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
\r
1543 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
\r
1545 else if (setVer==1) {
\r
1546 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1547 AddLayer((char*)"its2", 2.8 , x0, xrho, resRPhi, resZ);
\r
1548 AddLayer((char*)"its3", 3.6 , x0, xrho, resRPhi, resZ);
\r
1549 AddLayer((char*)"its4", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1550 AddLayer((char*)"its5", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1551 AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1552 AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1556 KMCTrackSummary::Bits(1,1,1),
\r
1557 KMCTrackSummary::Bits(0,0,0,1,1),
\r
1558 KMCTrackSummary::Bits(0,0,0,0,0,1,1)
\r
1560 RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
\r
1563 else if (setVer==2) {
\r
1564 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1565 AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ);
\r
1566 AddLayer((char*)"its2", 3.6 , x0, xrho, resRPhi, resZ);
\r
1567 AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ);
\r
1568 AddLayer((char*)"its3", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1569 AddLayer((char*)"its4", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1570 AddLayer((char*)"its5", 33.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1571 AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1572 AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ);
\r
1576 KMCTrackSummary::Bits(1,1,1,1),
\r
1577 KMCTrackSummary::Bits(0,0,0,0,1,1),
\r
1578 KMCTrackSummary::Bits(0,0,0,0,0,0,1,1,1)
\r
1580 RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
\r
1584 // last 2 layers strips
\r
1585 double resRPStr=0.0020, resZStr = 0.0830;
\r
1586 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1587 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1588 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1589 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1590 AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
\r
1591 AddLayer((char*)"its6", 39.6 , x0, xrho, resRPStr, resZStr);
\r
1592 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPStr, resZStr);
\r
1596 // resolution scaled with R as res(2.2) * R/2.2
\r
1597 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1598 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1599 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1600 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1601 AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
\r
1602 AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
\r
1603 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
\r
1604 KMCLayer* lr0 = 0;
\r
1605 for (int i=0;i<GetNActiveITSLayers();i++) {
\r
1606 KMCLayer *lr = GetActiveLayer(i);
\r
1607 if (lr->IsVertex()) continue;
\r
1612 double scl = 5*lr->GetRadius()/lr0->GetRadius();
\r
1613 SetResolution((char*)lr->GetName(), resRPhi*scl, resZ*scl*4);
\r
1617 // 1st 2 layers double
\r
1618 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1619 AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ);
\r
1621 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1622 AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ);
\r
1624 AddLayer((char*)"its3a", 6.4 , x0, xrho, resRPhi, resZ);
\r
1625 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1627 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1628 AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
\r
1629 AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
\r
1630 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
\r
1633 // last 4 layers doubled
\r
1634 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1635 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1636 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1638 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1639 AddLayer((char*)"its4a", 12.8 , x0, xrho, resRPhi, resZ);
\r
1641 AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
\r
1642 AddLayer((char*)"its5a", 23.9 , x0, xrho, resRPhi, resZ);
\r
1644 AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
\r
1645 AddLayer((char*)"its6a", 40.0 , x0, xrho, resRPhi, resZ);
\r
1647 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
\r
1648 AddLayer((char*)"its7a", 43.4 , x0, xrho, resRPhi, resZ);
\r
1651 /* //last 3 lr together
\r
1652 AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
\r
1653 AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
\r
1654 AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
\r
1655 AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
\r
1656 AddLayer((char*)"its5", 42.2 , x0, xrho, resRPhi, resZ);
\r
1657 AddLayer((char*)"its6", 42.6 , x0, xrho, resRPhi, resZ);
\r
1658 AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
\r
1661 AddTPC(0.1,0.1); // TPC
\r
1666 void KMCDetector::MakeAliceCurrent(Bool_t flagTPC, Int_t AlignResiduals) {
\r
1668 // Numbers taken from
\r
1669 // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
\r
1670 // number for misalingment: private communication with Andrea Dainese
\r
1672 // 1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};
\r
1674 AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation
\r
1675 if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
\r
1676 printf("vertex will work as constraint: %.4f %.4f",fgVtxConstraint[0],fgVtxConstraint[1]);
\r
1677 SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
\r
1680 AddLayer((char*)"bpipe_its",2.94,0.0022, 1.48e-01); // beam pipe
\r
1681 AddLayer((char*)"tshld1_its",11.5,0.0065, 1.34e-01); // Thermal shield // 1.3% /2
\r
1682 AddLayer((char*)"tshld2_its",31.0,0.0108, 2.22e-01); // Thermal shield // 1.3% /2
\r
1685 AddTPC(0.1,0.1); // TPC
\r
1687 // Adding the ITS - current configuration
\r
1689 if (AlignResiduals==0) {
\r
1690 AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, 0.0012, 0.0130);
\r
1691 AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, 0.0012, 0.0130);
\r
1692 AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, 0.0035, 0.0025);
\r
1693 AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, 0.0035, 0.0025);
\r
1694 AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 0.0020, 0.0830);
\r
1695 AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 0.0020, 0.0830);
\r
1698 KMCTrackSummary::Bits(1,1),
\r
1699 KMCTrackSummary::Bits(0,0,1,1),
\r
1700 KMCTrackSummary::Bits(0,0,0,0,1,1)
\r
1702 RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
\r
1704 } else if (AlignResiduals==1) {
\r
1706 // tracking errors ...
\r
1707 // (Additional systematic errors due to misalignments) ...
\r
1708 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
\r
1709 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
\r
1711 AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
\r
1712 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1713 AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
\r
1714 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1715 AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
\r
1716 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1717 AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
\r
1718 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1719 AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
\r
1720 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1721 AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
\r
1722 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1724 } else if (AlignResiduals==2) {
\r
1727 // tracking errors ... PLUS ... module misalignment
\r
1729 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
\r
1730 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
\r
1732 // the ITS modules are misalignment with small gaussian smearings with
\r
1733 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
\r
1735 AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
\r
1736 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1737 AddLayer((char*)"spd2_ots", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
\r
1738 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1739 AddLayer((char*)"sdd1_ots",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
\r
1740 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1741 AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
\r
1742 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1743 AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
\r
1744 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1745 AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
\r
1746 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1750 // the ITS modules are misalignment with small gaussian smearings with
\r
1751 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
\r
1752 // unknown in Z ????
\r
1754 AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
\r
1755 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
\r
1756 AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
\r
1757 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
\r
1758 AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
\r
1759 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
\r
1760 AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
\r
1761 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
\r
1762 AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
\r
1763 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
\r
1764 AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
\r
1765 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
\r
1771 void KMCDetector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
\r
1773 // Produces the standard performace plots
\r
1778 TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);
\r
1781 c1->cd(1); gPad->SetGridx(); gPad->SetGridy();
\r
1783 TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);
\r
1784 eff->SetTitle("Efficiencies");
\r
1786 if (!onlyPionEff) {
\r
1787 GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
\r
1788 GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
\r
1790 c1->cd(2); gPad->SetGridx(); gPad->SetGridy();
\r
1791 gPad->SetLogy(); gPad->SetLogx();
\r
1792 GetGraphMomentumResolution(color,linewidth)->Draw("AL");
\r
1794 c1->cd(3); gPad->SetGridx(); gPad->SetGridy();
\r
1796 GetGraphPointingResolution(0,color,linewidth)->Draw("AL");
\r
1798 c1->cd(4); gPad->SetGridx(); gPad->SetGridy();
\r
1800 GetGraphPointingResolution(1,color,linewidth)->Draw("AL");
\r
1804 TVirtualPad *c1 = gPad->GetMother();
\r
1807 GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");
\r
1808 if (!onlyPionEff) {
\r
1809 GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
\r
1810 GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
\r
1812 c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");
\r
1814 c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");
\r
1816 c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");
\r
1822 void KMCDetector::ApplyMS(KMCProbe* trc, double x2X0) const
\r
1824 // simulate random modification of track params due to the MS
\r
1825 if (x2X0<=0) return;
\r
1826 double alpha = trc->GetAlpha(); // store original alpha
\r
1827 double mass = trc->GetMass();
\r
1829 double snp = trc->GetSnp();
\r
1830 double dip = trc->GetTgl();
\r
1831 Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));
\r
1834 static double covCorr[15],covDum[21]={0};
\r
1835 static double mom[3],pos[3];
\r
1836 double *cov = (double*) trc->GetCovariance();
\r
1837 memcpy(covCorr,cov,15*sizeof(double));
\r
1839 trc->GetPxPyPz(mom);
\r
1840 double pt2 = mom[0]*mom[0]+mom[1]*mom[1];
\r
1841 double pt = TMath::Sqrt(pt2);
\r
1842 double ptot2 = pt2 + mom[2]*mom[2];
\r
1843 double ptot = TMath::Sqrt(ptot2);
\r
1844 double beta = ptot/TMath::Sqrt(ptot2 + mass*mass);
\r
1845 double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);
\r
1848 double phiSC = gRandom->Rndm()*TMath::Pi();
\r
1849 double thtSC = gRandom->Gaus(0,1.4142*sigth);
\r
1850 // printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);
\r
1851 double sn = TMath::Sin(thtSC);
\r
1852 double dx = sn*TMath::Sin(phiSC);
\r
1853 double dy = sn*TMath::Cos(phiSC);
\r
1854 double dz = TMath::Cos(thtSC);
\r
1856 // printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);
\r
1857 for (int i=3;i--;) mom[i] /= ptot;
\r
1858 double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
\r
1859 if (!IsZero(pt)) {
\r
1860 double pd1 = mom[0]/vmm;
\r
1861 double pd2 = mom[1]/vmm;
\r
1862 v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;
\r
1863 v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;
\r
1864 v[2] = -vmm*dx + mom[2]*dz;
\r
1869 v[2] = dz*TMath::Sign(1.,mom[2]);
\r
1871 double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
\r
1872 // printf("before :%+e %+e %+e || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2], sigth, x2X0, pt, beta);
\r
1874 // direction cosines -> p
\r
1875 for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;
\r
1876 // printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);
\r
1877 trc->Set(pos,mom,covDum,trc->Charge());
\r
1879 trc->Rotate(alpha);
\r
1880 memcpy(cov,covCorr,15*sizeof(double));
\r
1884 //________________________________________________________________________________
\r
1885 Bool_t KMCDetector::TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const
\r
1887 // Transport track till layer maxLr, applying random MS
\r
1889 for (Int_t j=0; j<maxLr; j++) {
\r
1890 KMCLayer* lr0 = (KMCLayer*)fLayers.At(j);
\r
1891 KMCLayer* lr = (KMCLayer*)fLayers.At(j+1);
\r
1892 if (!lr0->IsVertex()) {
\r
1893 ApplyMS(probTr,lr0->fx2X0); // apply MS
\r
1894 if (!probTr->CorrectForMeanMaterial(lr0,kFALSE)) return kFALSE;
\r
1897 if (!PropagateToLayer(probTr,lr,1)) return kFALSE;
\r
1898 // store randomized cluster local coordinates and phi
\r
1899 lr->ResetBgClusters();
\r
1901 gRandom->Rannor(rz,ry);
\r
1902 lr->GetMCCluster()->Set(probTr->GetY()+ry*lr->GetPhiRes(),probTr->GetZ()+rz*lr->GetZRes(),
\r
1903 probTr->GetX(), probTr->GetAlpha() );
\r
1910 //________________________________________________________________________________
\r
1911 Bool_t KMCDetector::SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr,int nMC, int offset)
\r
1913 // analityc and fullMC (nMC trials) evaluaion of tracks with given kinematics.
\r
1914 // the results are filled in KMCTrackSummary objects provided via summArr array
\r
1916 int progressP = 10;//0; // print progress in percents
\r
1918 progressP = int(nMC*0.01*progressP);
\r
1920 if (!SolveSingleTrackViaKalman(mass,pt,eta)) return kFALSE;
\r
1922 // Store non-updated track errors of inward propagated seed >>>>>>>>
\r
1923 int maxLr = fLastActiveITSLayer + offset;
\r
1924 if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
\r
1925 KMCProbe probeTmp = fProbe; // original probe at vertex
\r
1927 for (Int_t j=1; j<=maxLr; j++) {
\r
1929 // printf("Here0: %d\n",j);
\r
1930 if (!PropagateToLayer(&probeTmp,lr,1)) return 0;
\r
1931 if (j!=maxLr) if (!probeTmp.CorrectForMeanMaterial(lr, kFALSE)) return 0;
\r
1932 // printf("Prelim. Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(probeTmp.GetSigmaY2()),TMath::Sqrt(probeTmp.GetSigmaZ2()));
\r
1934 for (Int_t j=maxLr; j>0; j--) {
\r
1936 // printf("Here1: %d\n",j);
\r
1937 if (j!=maxLr) if (!PropagateToLayer(&probeTmp,lr,-1)) return 0;
\r
1938 lr->fSig2EstD = probeTmp.GetSigmaY2();
\r
1939 lr->fSig2EstZ = probeTmp.GetSigmaZ2();
\r
1940 // probeTmp.Print("l");
\r
1941 printf("Natural Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(lr->fSig2EstD),TMath::Sqrt(lr->fSig2EstZ));
\r
1942 if (!probeTmp.CorrectForMeanMaterial(lr, kTRUE)) return 0;
\r
1944 // Store non-updated track errors of inward propagated seed <<<<<<<<
\r
1946 int nsm = sumArr ? sumArr->GetEntriesFast() : 0;
\r
1947 KMCLayer* vtx = GetLayer(0);
\r
1949 for (int i=0;i<nsm;i++) {
\r
1950 KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(i);
\r
1951 if (!tsm) continue;
\r
1952 tsm->SetRefProbe( GetProbeTrack() ); // attach reference track (generated)
\r
1953 tsm->SetAnProbe( vtx->GetAnProbe() ); // attach analitycal solution
\r
1958 for (int it=0;it<nMC;it++) {
\r
1959 printf("ev: %d\n",it);
\r
1960 SolveSingleTrackViaKalmanMC(offset);
\r
1961 KMCProbe* trc = vtx->GetWinnerMCTrack();
\r
1962 vtx->GetMCTracks()->Print();
\r
1963 if (progressP==1 || (progressP>0 && (it%progressP)==0)) {
\r
1964 printf("%d%% done |",it*100/nMC);
\r
1965 sw.Stop(); sw.Print(); sw.Start(kFALSE);
\r
1967 for (int ism=nsm;ism--;) { // account the track in each of summaries
\r
1968 KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(ism);
\r
1969 if (!tsm) continue;
\r
1970 tsm->AddUpdCalls(GetUpdCalls());
\r
1971 tsm->AddTrack(trc);
\r
1976 printf("Total time: "); sw.Print();
\r
1980 //________________________________________________________________________________
\r
1981 Int_t KMCDetector::GetLayerID(int actID) const
\r
1983 // find physical layer id from active id
\r
1984 if (actID<0 || actID>fNActiveLayers) return -1;
\r
1985 int start = actID<fNActiveITSLayers ? fLastActiveITSLayer : fLastActiveLayer;
\r
1986 for (int i=start+1;i--;) {
\r
1987 if (GetLayer(i)->GetActiveID()==actID) return i;
\r
1992 //________________________________________________________________________________
\r
1993 KMCProbe* KMCDetector::KalmanSmooth(int actLr, int actMin,int actMax) const
\r
1995 // estimate kalman smoothed track params at given active lr
\r
1996 // from fit at layers actMin:actMax (excluding actLr)
\r
1997 // SolveSingleTrackViaKalman must have been called before
\r
1999 if (actMin>actMax) swap(actMin,actMax);
\r
2000 if (actMax>=fNActiveLayers) actMax = fNActiveLayers-1;
\r
2001 int nlrfit = actMax-actMin;
\r
2002 if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
\r
2003 if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
\r
2004 static KMCProbe iwd,owd;
\r
2006 // find phisical layer id's
\r
2007 int pLr = GetLayerID(actLr);
\r
2008 int pMin = GetLayerID(actMin);
\r
2009 int pMax = GetLayerID(actMax);
\r
2011 // printf(">>> %d %d %d\n",pLr, pMin,pMax);
\r
2012 Bool_t useIwd=kFALSE, useOwd=kFALSE;
\r
2013 if (pLr<pMax) { // need inward piece
\r
2014 iwd = GetLayer(pMax)->fTrCorr;
\r
2015 iwd.ResetCovMat();
\r
2016 iwd.GetHitsPatt() = 0;
\r
2017 for (int i=pMax;i>=pLr;i--) {
\r
2018 KMCLayer* lr = GetLayer(i);
\r
2019 // printf("IWD %d\n",i);
\r
2020 if (!lr->IsDead() && i!=pLr && i>=pMin) if (!UpdateTrack(&iwd,lr,&lr->fClCorr)) return 0;
\r
2022 if (!iwd.CorrectForMeanMaterial(lr,kTRUE)) return 0; // correct for materials of this layer
\r
2023 if (!PropagateToLayer(&iwd,GetLayer(i-1),-1)) return 0; // propagate to next layer
\r
2025 // printf("IWD%d: ",i); iwd.Print("l");
\r
2029 if (pLr>pMin) { // need outward piece
\r
2030 owd = GetLayer(pMin)->fTrCorr;
\r
2031 owd.ResetCovMat();
\r
2032 owd.GetHitsPatt() = 0;
\r
2033 for (int i=pMin;i<=pLr;i++) {
\r
2034 KMCLayer* lr = GetLayer(i);
\r
2035 // printf("OWD %d\n",i);
\r
2036 if (!lr->IsDead() && i!=pLr && i<=pMax) if (!UpdateTrack(&owd,lr,&lr->fClCorr)) return 0;
\r
2038 if (!owd.CorrectForMeanMaterial(lr,0)) return 0; // correct for materials of this layer
\r
2039 if (!PropagateToLayer(&owd,GetLayer(i+1), 1)) return 0; // propagate to next layer
\r
2041 // printf("OWD%d: ",i); owd.Print("l");
\r
2046 // was this extrapolation outside the fit range?
\r
2047 if (!useIwd) return (KMCProbe*)&owd;
\r
2048 if (!useOwd) return (KMCProbe*)&iwd;
\r
2050 // weight both tracks
\r
2051 if (!iwd.Propagate(owd.GetAlpha(),owd.GetX(),fBFieldG)) return 0;
\r
2052 double meas[2] = {owd.GetY(),owd.GetZ()};
\r
2053 double measErr2[3] = {owd.GetSigmaY2(), owd.GetSigmaZY(), owd.GetSigmaZ2()};
\r
2054 // printf("Weighting\n");
\r
2055 // owd.Print("l");
\r
2056 // iwd.Print("l");
\r
2057 if (!iwd.Update(meas,measErr2)) return 0;
\r
2058 iwd.GetHitsPatt() |= owd.GetHitsPatt();
\r
2060 // printf("->\n");
\r
2061 // iwd.Print("l");
\r
2063 return (KMCProbe*)&iwd;
\r
2067 //________________________________________________________________________________
\r
2068 KMCProbe* KMCDetector::KalmanSmoothFull(int actLr, int actMin,int actMax) const
\r
2070 // estimate kalman smoothed track params at given active lr
\r
2071 // from fit at layers actMin:actMax (excluding actLr)
\r
2072 // SolveSingleTrackViaKalman must have been called before
\r
2074 static TClonesArray prediction("KMCProbe",10);
\r
2075 static TClonesArray update("KMCProbe",10);
\r
2076 static KMCProbe res;
\r
2078 if (actMin>actMax) swap(actMin,actMax);
\r
2079 int nlrfit = actMax-actMin;
\r
2080 if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
\r
2081 if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
\r
2083 // find phisical layer id's
\r
2084 int pLr = GetLayerID(actLr);
\r
2085 int pMin = GetLayerID(actMin);
\r
2086 int pMax = GetLayerID(actMax);
\r
2088 int dir=0,dirInt=0;
\r
2089 if (pLr<=pMin) dir=-1; // inward extrapolation
\r
2090 else if (pLr>=pMax) dir= 1; // outward extrapolation
\r
2091 else if (actMax-actLr >= actLr-actMin) dirInt = -1; // inward interpolation (the test point is closer to inner layer)
\r
2092 else dirInt = 1; // outward interpolation (the test point is closer to outer layer)
\r
2094 if (dir!=0) { // no sens to do smoothing: simple Kalman filtering extrapolation
\r
2095 int start = dir<0 ? pMax : pMin;
\r
2096 res = GetLayer(start)->fTrCorr;
\r
2097 res.ResetCovMat();
\r
2099 for (int i=(dir<0?pMax:pMin); i!=pLr; i+=dir) { // track till nearest layer to pLr
\r
2101 if (!lr->IsDead() && !(i<pMin ||i>pMax)) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range
\r
2102 if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer
\r
2103 if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer
\r
2105 if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this nearest layer
\r
2106 if (!PropagateToLayer(&res,GetLayer(pLr), dir)) return 0; // propagate to test layer
\r
2107 return (KMCProbe*)&res;
\r
2110 // too bad, need to do real filtering
\r
2112 int start = dirInt<0 ? pMax : pMin;
\r
2113 int stop = dirInt<0 ? pMin-1 : pMax+1;
\r
2114 res = GetLayer(start)->fTrCorr;
\r
2115 res.ResetCovMat();
\r
2118 for (int i=start; i!=stop; i+=dirInt) { // track in full range, storing updates and predictions
\r
2119 new(prediction[count]) KMCProbe(res);
\r
2121 if (!lr->IsDead() && i!=pLr) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range
\r
2122 new(update[count]) KMCProbe(res);
\r
2123 if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer
\r
2124 if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer
\r
2127 return (KMCProbe*)&res;
\r
2131 //________________________________________________________________________________
\r
2132 Bool_t KMCDetector::SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta)
\r
2134 // analytical estimate of tracking resolutions
\r
2135 // fProbe.SetUseLogTermMS(kTRUE);
\r
2137 if (fMinITSHits>fNActiveITSLayers) {fMinITSHits = fNActiveITSLayers; printf("Redefined request of min N ITS hits to %d\n",fMinITSHits);}
\r
2138 if (TMath::Abs(eta)<1e-3) fDensFactorEta = 1.;
\r
2140 fDensFactorEta = TMath::Tan( 2.*TMath::ATan(TMath::Exp(-TMath::Abs(eta))) );
\r
2141 fDensFactorEta = 1./TMath::Sqrt( 1. + 1./fDensFactorEta/fDensFactorEta);
\r
2143 double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta));
\r
2144 KMCProbe* probe = PrepareKalmanTrack(pt,lambda,mass,-1);
\r
2145 if (!probe) return kFALSE;
\r
2150 // Start the track fitting --------------------------------------------------------
\r
2152 // Back-propagate the covariance matrix along the track.
\r
2153 // Kalman loop over the layers
\r
2155 KMCProbe* currTr = 0;
\r
2156 lr = (KMCLayer*)fLayers.At(fLastActiveLayerTracked);
\r
2157 lr->fTrCorr = *probe;
\r
2158 delete probe; // rethink...
\r
2160 for (Int_t j=fLastActiveLayerTracked; j--; ) { // Layer loop
\r
2162 KMCLayer *lrP = lr;
\r
2163 lr = (KMCLayer*)fLayers.At(j);
\r
2165 lr->fTrCorr = lrP->fTrCorr;
\r
2166 currTr = &lr->fTrCorr;
\r
2167 currTr->ResetHit(lrP->GetActiveID());
\r
2169 // if there was a measurement on prev layer, update the track
\r
2170 if (!lrP->IsDead()) { // include measurement
\r
2171 KMCCluster cl(currTr->GetY(),currTr->GetZ(), currTr->GetX(), currTr->GetAlpha());
\r
2172 if (!UpdateTrack(currTr,lrP,&cl)) return kFALSE;
\r
2174 if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) return kFALSE; // correct for materials of this layer
\r
2175 if (!PropagateToLayer(currTr,lr,-1)) return kFALSE; // propagate to current layer
\r
2177 } // end loop over layers
\r
2182 //____________________________________________________________
\r
2183 Bool_t KMCDetector::SolveSingleTrackViaKalmanMC(int offset)
\r
2185 // MC estimate of tracking resolutions/effiencies. Requires that the SolveSingleTrackViaKalman
\r
2186 // was called before, since it uses data filled by this method
\r
2188 // The MC tracking will be done starting from fLastActiveITSLayer + offset (before analytical estimate will be used)
\r
2190 // At this point, the fProbe contains the track params generated at vertex.
\r
2191 // Clone it and propagate to target layer to generate hit positions affected by MS
\r
2194 KMCProbe *currTrP=0,*currTr=0;
\r
2195 int maxLr = fLastActiveITSLayer + offset;
\r
2196 if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
\r
2197 ResetMCTracks(maxLr);
\r
2198 KMCLayer* lr = (KMCLayer*)fLayers.At(maxLr);
\r
2199 currTr = lr->AddMCTrack(&fProbe); // start with original track at vertex
\r
2201 if (!TransportKalmanTrackWithMS(currTr, maxLr)) return kFALSE; // transport it to outermost layer where full MC is done
\r
2203 if (fLastActiveITSLayer<fLastActiveLayerTracked) { // prolongation from TPC
\r
2204 // start from correct track propagated from above till maxLr
\r
2205 double *covMS = (double*)currTr->GetCovariance();
\r
2206 const double *covIdeal =lr->fTrCorr.GetCovariance();
\r
2207 for (int i=15;i--;) covMS[i] = covIdeal[i];
\r
2209 else { // ITS SA: randomize the starting point
\r
2210 // double *pars = (double*)currTr->GetParameter();
\r
2211 // pars[0] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaY2()));
\r
2212 // pars[1] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaZ2()));
\r
2214 currTr->ResetCovMat();
\r
2216 double *trCov = (double*)currTr->GetCovariance();
\r
2217 double *trPars = (double*)currTr->GetParameter();
\r
2218 const double kLargeErr2PtI = 0.3*0.3;
\r
2219 trCov[14] = TMath::Max(trCov[14],kLargeErr2PtI*trPars[4]*trPars[4]);
\r
2223 for (Int_t j=maxLr; j--; ) { // Layer loop
\r
2225 KMCLayer *lrP = lr;
\r
2226 lr = (KMCLayer*)fLayers.At(j);
\r
2227 int ntPrev = lrP->GetNMCTracks();
\r
2229 if (lrP->IsDead()) { // for passive layer just propagate the copy of all tracks of prev layer >>>
\r
2230 for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
\r
2231 currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
\r
2232 currTr = lr->AddMCTrack( currTrP );
\r
2233 if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
\r
2234 if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer
\r
2237 } // treatment of dead layer <<<
\r
2239 if (lrP->IsTPC()) { // we don't consider bg hits in TPC, just update with MC cluster
\r
2240 for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
\r
2241 currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
\r
2242 currTr = lr->AddMCTrack( currTrP );
\r
2243 if (!UpdateTrack(currTr, lrP, lrP->GetMCCluster(), kTRUE)) {currTr->Kill(); continue;} // update with correct MC cl.
\r
2244 if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
\r
2245 if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer
\r
2248 } // treatment of ideal (TPC?) layer <<<
\r
2250 // active layer under eff. study (ITS?): propagate copy of every track to MC cluster frame (to have them all in the same frame)
\r
2251 // and calculate the limits of bg generation
\r
2252 KMCCluster* clMC = lrP->GetMCCluster();
\r
2253 if (lrP->GetLayerEff()<gRandom->Rndm()) clMC->Kill(); // simulate inefficiency
\r
2254 ResetSearchLimits();
\r
2256 for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
\r
2257 currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
\r
2258 currTr = lr->AddMCTrack( currTrP );
\r
2259 if (!currTr->PropagateToCluster(clMC,fBFieldG)) {currTr->Kill(); continue;} // go to MC cluster
\r
2260 if ( !(currTr->GetNITSHits()>0 && currTr->GetNITSHits()==currTr->GetNFakeITSHits()) ) UpdateSearchLimits(currTr, lrP); // RS
\r
2264 // printf("%3d seeds\n",nseeds);
\r
2265 if (fUseBackground && lrP->IsITS()) GenBgClusters(lrP); // generate background hits
\r
2267 ntPrev = lr->GetNMCTracks();
\r
2268 for (int itr=ntPrev;itr--;) { // loop over all tracks PROPAGATED from previous layer to clusters frame on previous layer
\r
2269 currTrP = lr->GetMCTrack(itr); // this is a seed from prev layer. The new clusters are attached to its copies, the seed itself
\r
2270 // will be propagated w/o cluster update if it does not violate requested "reconstructed" track settings
\r
2271 if (currTrP->IsKilled()) continue;
\r
2272 //printf("Check %d %p %d\n",itr,currTrP,currTrP->GetUniqueID()); currTrP->Print();
\r
2273 CheckTrackProlongations(currTrP, lr, lrP);
\r
2274 if (NeedToKill(currTrP)) currTrP->Kill(); // kill track which was not updated at lrP
\r
2275 //currTrP->Kill(); // kill track which was not updated at lrP
\r
2278 lr->GetMCTracks()->Sort();
\r
2279 int ntTot = lr->GetNMCTracks(); // propagate max amount of allowed tracks to current layer
\r
2280 if (ntTot>fMaxSeedToPropagate && fMaxSeedToPropagate>0) {
\r
2281 for (int itr=ntTot;itr>=fMaxSeedToPropagate;itr--) lr->GetMCTracks()->RemoveAt(itr);
\r
2282 ntTot = fMaxSeedToPropagate;
\r
2285 for (int itr=ntTot;itr--;) {
\r
2286 currTr = lr->GetMCTrack(itr);
\r
2287 if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill();continue;} // correct for materials of prev. layer
\r
2288 if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill();continue;} // propagate to current layer
\r
2290 AliDebug(1,Form("Got %d tracks on layer %s",ntTot,lr->GetName()));
\r
2291 // lr->GetMCTracks()->Print();
\r
2293 } // end loop over layers
\r
2295 // do we use vertex constraint?
\r
2296 KMCLayer *vtx = GetLayer(0);
\r
2297 if (!vtx->IsDead() && vtx->IsITS()) {
\r
2298 int ntr = vtx->GetNMCTracks();
\r
2299 for (int itr=0;itr<ntr;itr++) {
\r
2300 currTr = vtx->GetMCTrack(itr);
\r
2301 if (currTr->IsKilled()) continue;
\r
2302 KMCCluster* clv = vtx->GetMCCluster();
\r
2303 double meas[2] = {clv->GetY(),clv->GetZ()};
\r
2304 double measErr2[3] = {vtx->fPhiRes*vtx->fPhiRes,0,vtx->fZRes*vtx->fZRes};
\r
2305 double chi2v = currTr->GetPredictedChi2(meas,measErr2);
\r
2306 currTr->AddHit(vtx->GetActiveID(), chi2v, -1);
\r
2307 currTr->SetInnerLrChecked(vtx->GetActiveID());
\r
2308 if (NeedToKill(currTr)) currTr->Kill();
\r
2309 // if (vtx->IsITS()) {if (!UpdateTrack(currTr, vtx, vtx->GetMCCluster(), kFALSE)) {currTr->Kill();continue;}}
\r
2312 EliminateUnrelated();
\r
2317 //____________________________________________________________________________
\r
2318 Bool_t KMCDetector::PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const
\r
2320 // bring the track to layer and rotat to frame normal to its surface
\r
2321 if (!trc->PropagateToR(lr->fR,fBFieldG, dir)) return kFALSE;
\r
2323 // rotate to frame with X axis normal to the surface (defined by ideal track)
\r
2324 if (!lr->IsVertex()) {
\r
2326 trc->GetXYZ(pos); // lab position
\r
2327 double phi = TMath::ATan2(pos[1],pos[0]);
\r
2328 if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;
\r
2329 if (!trc->Rotate(phi)) {
\r
2330 AliDebug(2,Form("Failed to rotate to the frame (phi:%+.3f) of layer at %.2f at XYZ: %+.3f %+.3f %+.3f",
\r
2331 phi,lr->fR,pos[0],pos[1],pos[2]));
\r
2332 if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
\r
2340 //____________________________________________________________________________
\r
2341 Bool_t KMCDetector::UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster) const
\r
2343 // update track with measured cluster
\r
2344 // propagate to cluster
\r
2345 double meas[2] = {cl->GetY(),cl->GetZ()}; // ideal cluster coordinate
\r
2346 double measErr2[3] = {lr->fPhiRes*lr->fPhiRes,0,lr->fZRes*lr->fZRes};
\r
2348 if (goToCluster) if (!trc->PropagateToCluster(cl,fBFieldG)) return kFALSE; // track was not propagated to cluster frame
\r
2350 double chi2 = trc->GetPredictedChi2(meas,measErr2);
\r
2351 // if (chi2>fMaxChi2Cl) return kTRUE; // chi2 is too large
\r
2353 if (!trc->Update(meas,measErr2)) {
\r
2354 AliDebug(2,Form("layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
\r
2355 lr->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
\r
2356 if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
\r
2359 trc->AddHit(lr->GetActiveID(), chi2);
\r
2364 //____________________________________________________________________________
\r
2365 Int_t KMCDetector::GenBgClusters(KMCLayer* lr)
\r
2367 // Generate fake clusters in precalculated RPhi,Z range
\r
2368 if (fNBgLimits<1) return 0; // limits were not set - no track was prolongated
\r
2370 // Fix search limits to avoid seeds which will anyway point very far from the vertex
\r
2371 double tolY = TMath::Sqrt(lr->fSig2EstD)*fMaxChi2ClSQ;
\r
2372 double tolZ = TMath::Sqrt(lr->fSig2EstZ)*fMaxChi2ClSQ;
\r
2374 // printf("Before: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
\r
2375 if (fBgYMin < lr->fClCorr.fY-tolY) fBgYMin = lr->fClCorr.fY-tolY;
\r
2376 if (fBgYMax > lr->fClCorr.fY+tolY) fBgYMax = lr->fClCorr.fY+tolY;
\r
2377 if (fBgZMin < lr->fClCorr.fZ-tolZ) fBgZMin = lr->fClCorr.fZ-tolZ;
\r
2378 if (fBgZMax > lr->fClCorr.fZ+tolZ) fBgZMax = lr->fClCorr.fZ+tolZ;
\r
2379 //printf("After: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
\r
2381 double dy = fBgYMax - fBgYMin;
\r
2382 double dz = fBgZMax - fBgZMin;
\r
2383 double surf = dy*dz; // surface of generation
\r
2384 if (surf<0) return 0;
\r
2385 double poissProb = surf*HitDensity(lr->fR)*lr->GetLayerEff();
\r
2386 AliDebug(2,Form("Bg for Lr %s (r=%.2f) : Density %.2f on surface %.2e [%+.4f : %+.4f][%+.4f %+.4f]",
\r
2387 lr->GetName(),lr->fR,HitDensity(lr->fR),surf,fBgYMin,fBgYMax,fBgZMin,fBgZMax));
\r
2388 int nFakesGen = gRandom->Poisson( poissProb ); // preliminary number of extra clusters to test
\r
2389 KMCCluster *refCl = lr->GetMCCluster();
\r
2390 double sig2y = lr->GetPhiRes()*lr->GetPhiRes();
\r
2391 double sig2z = lr->GetZRes()*lr->GetZRes();
\r
2392 for (int ic=nFakesGen;ic--;) {
\r
2393 double y = fBgYMin+dy*gRandom->Rndm();
\r
2394 double z = fBgZMin+dz*gRandom->Rndm();
\r
2395 double dfy = y-refCl->GetY();
\r
2396 double dfz = z-refCl->GetZ();
\r
2397 double dist = (dfy*dfy)/sig2y + (dfz*dfz)/sig2z;
\r
2398 if (dist<4) continue; // avoid overlap with MC cluster
\r
2399 lr->AddBgCluster(y, z, refCl->GetX(), refCl->GetPhi());
\r
2401 AliDebug(2,Form("Added %6d noise clusters on lr %s (poisson Prob=%8.2f for surface %.2e) DY:%7.4f DZ: %7.4f",
\r
2402 lr->GetNBgClusters(),lr->GetName(),poissProb,surf,dy,dz));
\r
2407 //____________________________________________________________________________
\r
2408 void KMCDetector::UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr)
\r
2410 // define the search window for track on layer (where the bg hist will be generated)
\r
2411 static double *currYMin = fBgYMinTr.GetArray();
\r
2412 static double *currYMax = fBgYMaxTr.GetArray();
\r
2413 static double *currZMin = fBgZMinTr.GetArray();
\r
2414 static double *currZMax = fBgZMaxTr.GetArray();
\r
2416 double sizeY = probe->GetSigmaY2(), sizeZ = probe->GetSigmaZ2();
\r
2418 if (probe->GetNITSHits()<3) sizeY = 10*lr->fSig2EstD;
\r
2419 if (probe->GetNITSHits()<2) sizeZ = 10*lr->fSig2EstZ;
\r
2420 sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY+lr->fPhiRes*lr->fPhiRes); // max deviation in rphi to accept
\r
2421 sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ+lr->fZRes*lr->fZRes); // max deviation in dz to accept
\r
2425 if (probe->GetNITSHits()<3) sizeY = 1./(1./lr->fSig2EstD + 1./sizeY);
\r
2426 if (probe->GetNITSHits()<2) sizeZ = 1./(1./lr->fSig2EstZ + 1./sizeZ);
\r
2427 sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY); // max deviation in rphi to accept
\r
2428 sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ); // max deviation in dz to accept
\r
2431 // if (sizeY>2) sizeY=2;
\r
2432 // if (sizeZ>2) sizeZ=2;
\r
2433 // printf("Sizes at %s: %.5f %.5f\n",lr->GetName(), sizeY,sizeZ);
\r
2435 if (fNBgLimits>=fBgYMinTr.GetSize()) { // expand arrays, update pointers
\r
2436 fBgYMinTr.Set(2*(fNBgLimits+1));
\r
2437 fBgYMaxTr.Set(2*(fNBgLimits+1));
\r
2438 fBgZMinTr.Set(2*(fNBgLimits+1));
\r
2439 fBgZMaxTr.Set(2*(fNBgLimits+1));
\r
2440 currYMin = fBgYMinTr.GetArray();
\r
2441 currYMax = fBgYMaxTr.GetArray();
\r
2442 currZMin = fBgZMinTr.GetArray();
\r
2443 currZMax = fBgZMaxTr.GetArray();
\r
2445 if (fBgYMin > (currYMin[fNBgLimits]=probe->GetY()-sizeY) ) fBgYMin = currYMin[fNBgLimits];
\r
2446 if (fBgYMax < (currYMax[fNBgLimits]=probe->GetY()+sizeY) ) fBgYMax = currYMax[fNBgLimits];
\r
2447 if (fBgZMin > (currZMin[fNBgLimits]=probe->GetZ()-sizeZ) ) fBgZMin = currZMin[fNBgLimits];
\r
2448 if (fBgZMax < (currZMax[fNBgLimits]=probe->GetZ()+sizeZ) ) fBgZMax = currZMax[fNBgLimits];
\r
2449 if (AliLog::GetGlobalDebugLevel()>=2) {
\r
2450 probe->Print("l");
\r
2451 AliInfo(Form("Seed%3d Lr %s limits for y:%+8.4f z:%+8.4f [%+.4f : %+.4f][%+.4f %+.4f]",fNBgLimits,lr->GetName(),probe->GetY(),probe->GetZ(),currYMin[fNBgLimits],currYMax[fNBgLimits],currZMin[fNBgLimits],currZMax[fNBgLimits]));
\r
2452 AliInfo(Form("Global Limits Lr %s [%+.4f : %+.4f][%+.4f %+.4f]",lr->GetName(),fBgYMin,fBgYMax,fBgZMin,fBgZMax));
\r
2453 AliInfo(Form("MC Cluster: %+.4f : %+.4f",lr->fClMC.fY, lr->fClMC.fZ));
\r
2455 probe->SetUniqueID(fNBgLimits++);
\r
2457 if (lr->IsITS() && probe->GetNFakeITSHits()==0) {
\r
2458 if (fHMCLrResidRPhi) fHMCLrResidRPhi->Fill(probe->GetY() - lr->GetMCCluster()->GetY(), lr->GetActiveID());
\r
2459 if (fHMCLrResidZ) fHMCLrResidZ->Fill(probe->GetZ() - lr->GetMCCluster()->GetZ(),lr->GetActiveID());
\r
2464 //____________________________________________________________________________
\r
2465 void KMCDetector::CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP)
\r
2467 // explore prolongation of probe from lrP to lr with all possible clusters of lrP
\r
2468 // the probe is already brought to clusters frame
\r
2469 int nCl = lrP->GetNBgClusters();
\r
2470 double measErr2[3] = {lrP->fPhiRes*lrP->fPhiRes,0,lrP->fZRes*lrP->fZRes};
\r
2471 double meas[2] = {0,0};
\r
2472 UInt_t tmpID = probe->GetUniqueID();
\r
2473 double yMin = fBgYMinTr[tmpID];
\r
2474 double yMax = fBgYMaxTr[tmpID];
\r
2475 double zMin = fBgZMinTr[tmpID];
\r
2476 double zMax = fBgZMaxTr[tmpID];
\r
2478 probe->SetInnerLrChecked(lrP->GetActiveID());
\r
2479 for (int icl=-1;icl<nCl;icl++) {
\r
2480 KMCCluster* cl = icl<0 ? lrP->GetMCCluster() : lrP->GetBgCluster(icl); // -1 is for true MC cluster
\r
2481 if (cl->IsKilled()) {
\r
2482 if (AliLog::GetGlobalDebugLevel()>1) {printf("Skip cluster %d ",icl); cl->Print();}
\r
2485 double y = cl->GetY();
\r
2486 double z = cl->GetZ();
\r
2487 AliDebug(2,Form("Check seed%d against cl#%d out of %d at layer %s | y:%+8.4f z:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f]",tmpID,icl,nCl,lrP->GetName(),y,z,yMin,yMax,zMin,zMax));
\r
2488 if (AliLog::GetGlobalDebugLevel()>0) {
\r
2489 if (icl==-1 && probe->GetNFakeITSHits()==0) {
\r
2490 meas[0] = y; meas[1] = z;
\r
2491 double chi2a = probe->GetPredictedChi2(meas,measErr2);
\r
2492 if (chi2a>fMaxChi2Cl || (y<yMin || y>yMax) || (z<zMin || z>zMax)) {
\r
2494 printf("Loosing good point (y:%+8.4f z:%+8.4f) on lr %s: chi2: %.2f | dy:%+8.4f dz:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f] |x: %.2f %.2f | phi: %.2f %.2f\n",
\r
2495 y,z,lrP->GetName(),chi2a,y-probe->GetY(),z-probe->GetZ(),yMin,yMax,zMin,zMax, probe->GetX(), cl->GetX(), probe->GetAlpha(), cl->GetPhi());
\r
2499 if (y<yMin || y>yMax) continue; // preliminary check on Y
\r
2500 if (z<zMin || z>zMax) continue; // preliminary check on Z
\r
2501 meas[0] = y; meas[1] = z;
\r
2502 double chi2 = probe->GetPredictedChi2(meas,measErr2);
\r
2503 if (fHMCLrChi2 && probe->GetNFakeITSHits()==0 && icl==-1) fHMCLrChi2->Fill(chi2,lrP->GetActiveID());
\r
2504 AliDebug(2,Form("Seed-to-cluster chi2 = Chi2=%.2f",chi2));
\r
2505 if (chi2>fMaxChi2Cl) continue;
\r
2507 // update track copy
\r
2508 KMCProbe* newTr = lr->AddMCTrack( probe );
\r
2510 if (!newTr->Update(meas,measErr2)) {
\r
2511 AliDebug(2,Form("Layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
\r
2512 lrP->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
\r
2513 if (AliLog::GetGlobalDebugLevel()>1) newTr->Print("l");
\r
2517 newTr->AddHit(lrP->GetActiveID(), chi2, icl);
\r
2518 if (AliLog::GetGlobalDebugLevel()>1) {
\r
2519 AliInfo("Cloned updated track is:");
\r
2522 if (NeedToKill(newTr)) newTr->Kill();
\r
2527 //____________________________________________________________________________
\r
2528 void KMCDetector::ResetMCTracks(Int_t maxLr)
\r
2530 int nl = GetNLayers();
\r
2531 if (maxLr<0 || maxLr>=nl) maxLr = nl-1;
\r
2532 for (int i=maxLr+1;i--;) GetLayer(i)->ResetMCTracks();
\r
2535 //____________________________________________________________________________
\r
2536 Bool_t KMCDetector::NeedToKill(KMCProbe* probe) const
\r
2538 // check if the seed at given layer (last one where update was tried)
\r
2539 // still has chances to be reconstructed
\r
2540 const Bool_t kModeKillMiss = kFALSE;
\r
2542 Bool_t kill = kFALSE;
\r
2544 int il = probe->GetInnerLayerChecked();
\r
2545 int nITS = probe->GetNITSHits();
\r
2546 int nITSMax = nITS + il; // maximum it can have
\r
2547 if (nITSMax<fMinITSHits) {
\r
2550 } // has no chance to collect enough ITS hits
\r
2552 int ngr = fPattITS.GetSize();
\r
2553 if (ngr>0) { // check pattern
\r
2554 UInt_t patt = probe->GetHitsPatt();
\r
2555 // complete the layers not checked yet
\r
2556 for (int i=il;i--;) patt |= (0x1<<i);
\r
2557 for (int ig=ngr;ig--;)
\r
2558 if (!(((UInt_t)fPattITS[ig]) & patt)) {
\r
2565 if (nITS>2) { // check if smallest possible norm chi2/ndf is acceptable
\r
2566 double chi2min = probe->GetChi2();
\r
2567 if (kModeKillMiss) {
\r
2568 int nMiss = fNActiveITSLayers - probe->GetInnerLayerChecked() - nITS; // layers already missed
\r
2569 chi2min = nMiss*probe->GetMissingHitPenalty();
\r
2571 chi2min /= ((nITSMax<<1)-KMCProbe::kNDOF);
\r
2572 if (chi2min>fMaxNormChi2NDF) {
\r
2578 // loose vertex constraint
\r
2581 probe->GetZAt(0,fBFieldG,dst);
\r
2582 //printf("Zd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
\r
2583 if (TMath::Abs(dst)>10.) {
\r
2589 probe->GetYAt(0,fBFieldG,dst);
\r
2590 //printf("Dd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
\r
2591 if (TMath::Abs(dst)>10.) {
\r
2599 if (kill && AliLog::GetGlobalDebugLevel()>1 && probe->GetNFakeITSHits()==0) {
\r
2600 printf("Killing good seed, last upd layer was %d\n",probe->GetInnerLayerChecked());
\r
2601 probe->Print("l");
\r
2606 //_____________________________________________________________________
\r
2607 void KMCDetector::EliminateUnrelated()
\r
2609 // kill useless tracks
\r
2610 KMCLayer* lr = GetLayer(0);
\r
2611 int ntr = lr->GetNMCTracks();
\r
2613 for (int itr=0;itr<ntr;itr++) {
\r
2614 KMCProbe* probe = lr->GetMCTrack(itr);
\r
2615 if (probe->IsKilled()) continue;
\r
2616 if (probe->GetNITSHits()-probe->GetNFakeITSHits()<1) {probe->Kill(); continue;}
\r
2619 lr->GetMCTracks()->Sort();
\r
2620 const int kDump = 0;
\r
2622 printf("Valid %d out of %d\n",nval, ntr);
\r
2623 ntr = ntr>kDump ? kDump:0;
\r
2624 for (int itr=0;itr<ntr;itr++) {
\r
2625 lr->GetMCTrack(itr)->Print();
\r
2630 //_____________________________________________________________________
\r
2631 void KMCDetector::RequirePattern(UInt_t *patt, int groups)
\r
2633 if (groups<1) {fPattITS.Set(0); return;}
\r
2634 fPattITS.Set(groups);
\r
2635 for (int i=0;i<groups;i++) fPattITS[i] = patt[i];
\r
2638 //_____________________________________________________________________
\r
2639 void KMCDetector::CalcHardSearchLimits(double dzv)
\r
2643 zlims.Set(fNActiveITSLayers);
\r
2644 for (int il0=0;il0<fNActiveITSLayers;il0++) {
\r
2645 KMCLayer* lr0 = GetActiveLayer(il0);
\r
2646 double angZTol = dzv/lr0->GetRadius();
\r
2647 for (int il1=0;il1<fNActiveITSLayers;il1++) {
\r
2648 if (il1==il0) continue;
\r
2649 KMCLayer* lr1 = GetActiveLayer(il1);
\r
2650 double ztol = angZTol*TMath::Abs(lr0->GetRadius() - lr1->GetRadius());
\r
2651 if (ztol>zlims[il1]) zlims[il1] = ztol;
\r
2655 for (int il=0;il<fNActiveITSLayers;il++) printf("ZTol%d: %8.4f\n",il,zlims[il]);
\r
2658 //_______________________________________________________
\r
2659 double KMCDetector::PropagateBack(KMCProbe* trc)
\r
2661 static KMCProbe bwd;
\r
2663 bwd.ResetCovMat();
\r
2664 static double measErr2[3] = {0,0,0};
\r
2665 static double meas[2] = {0,0};
\r
2667 double chi2Tot = 0;
\r
2668 for (int il=1;il<=fLastActiveITSLayer;il++) {
\r
2669 KMCLayer* lr = GetLayer(il);
\r
2670 if (!PropagateToLayer(&bwd,lr,1)) return -1;
\r
2671 int aID = lr->GetActiveID();
\r
2672 if (aID>-1 && (icl=bwd.fClID[aID])>=-1) {
\r
2673 KMCCluster* clMC = icl<0 ? lr->GetMCCluster() : lr->GetBgCluster(icl);
\r
2674 if (!bwd.PropagateToCluster(clMC,fBFieldG)) return -1;
\r
2675 meas[0] = clMC->GetY(); meas[1] = clMC->GetZ();
\r
2676 measErr2[0] = lr->fPhiRes*lr->fPhiRes;
\r
2677 measErr2[2] = lr->fZRes*lr->fZRes;
\r
2678 double chi2a = bwd.GetPredictedChi2(meas,measErr2);
\r
2680 printf("Chis %d (cl%+3d): t2c: %6.3f tot: %6.3f\n",aID,icl,chi2a, chi2Tot);
\r
2681 bwd.Update(meas,measErr2);
\r
2682 bwd.AddHit(aID, chi2a, icl);
\r
2684 if (!bwd.CorrectForMeanMaterial(lr,kFALSE)) return -1;
\r
2690 //////////////////////////////////////////////////////////////////////////////////////////////
\r
2691 //--------------------------------------------------------------------------------------------
\r
2692 ClassImp(KMCTrackSummary)
\r
2694 Int_t KMCTrackSummary::fgSumCounter = 0;
\r
2696 //____________________________________________________________________________
\r
2697 KMCTrackSummary::KMCTrackSummary(const char* name,const char* title, int nlr) :
\r
2698 TNamed(name,name), fNITSLayers(nlr),
\r
2699 fPattITS(0),fPattITSFakeExcl(0),fPattITSCorr(0),
\r
2700 fHMCChi2(0),fHMCSigDCARPhi(0),fHMCSigDCAZ(0),fHMCSigPt(0),fHMCNCl(0),fHMCFakePatt(0),
\r
2701 fCountAcc(0),fCountTot(0),fUpdCalls(0),
\r
2702 fRefProbe(),fAnProbe()
\r
2704 // create summary structure for specific track conditions
\r
2707 SetMinMaxClITSFake();
\r
2708 SetMinMaxClITSCorr();
\r
2710 if (name==0 && title==0 && nlr==0) return; // default dummy contructor
\r
2712 enum {kNBinRes=1000};
\r
2713 const double kMaxChi2(10),kMaxResRPhi=0.1, kMaxResZ=0.1,kMaxResPt=0.3;
\r
2715 TString strN = name, strT = title;
\r
2716 if (strN.IsNull()) strN = Form("TrSum%d",fgSumCounter);
\r
2717 if (strT.IsNull()) strT = strN;
\r
2720 if (fNITSLayers<1) {
\r
2721 fNITSLayers=KMCProbe::GetNITSLayers();
\r
2722 if (fNITSLayers<1) {AliError("N ITS layer is not provided and not available from KMCProbe::GetNITSLayers()"); exit(1);}
\r
2723 AliInfo(Form("%s No N Layers provided, set %d",strN.Data(),fNITSLayers));
\r
2725 nlr = fNITSLayers;
\r
2729 nam = Form("%s_mc_chi2",strN.Data());
\r
2730 tit = Form("%s MC #chi^{2}",strT.Data());
\r
2731 fHMCChi2 = new TH1F(nam.Data(),tit.Data(),kNBinRes,0,kMaxChi2);
\r
2732 fHMCChi2->GetXaxis()->SetTitle("#chi^{2}");
\r
2733 fHMCChi2->Sumw2();
\r
2735 nam = Form("%s_mc_DCArphi",strN.Data());
\r
2736 tit = Form("%s MC #DeltaDCA r#phi",strT.Data());
\r
2737 fHMCSigDCARPhi = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResRPhi,kMaxResRPhi);
\r
2738 fHMCSigDCARPhi->GetXaxis()->SetTitle("#Delta r#phi");
\r
2739 fHMCSigDCARPhi->Sumw2();
\r
2741 nam = Form("%s_mc_DCAz",strN.Data());
\r
2742 tit = Form("%s MC #DeltaDCA Z",strT.Data());
\r
2743 fHMCSigDCAZ = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResZ,kMaxResZ);
\r
2744 fHMCSigDCAZ->GetXaxis()->SetTitle("#Delta Z");
\r
2745 fHMCSigDCAZ->Sumw2();
\r
2747 nam = Form("%s_mc_pt",strN.Data());
\r
2748 tit = Form("%s MC $Deltap_{T}/p_{T}",strT.Data());
\r
2749 fHMCSigPt = new TH1F(nam.Data(),tit.Data(),2*kNBinRes,-kMaxResPt,kMaxResPt);
\r
2750 fHMCSigPt->GetXaxis()->SetTitle("#Deltap_{T}/p_{T}");
\r
2751 fHMCSigPt->Sumw2();
\r
2753 nam = Form("%s_n_cl",strN.Data());
\r
2754 tit = Form("%s N Clusters",strT.Data());
\r
2755 fHMCNCl = new TH2F(nam.Data(),tit.Data(),nlr,1,nlr+1,nlr,1,nlr+1);
\r
2756 fHMCNCl->GetXaxis()->SetTitle("N Clusters Tot");
\r
2757 fHMCNCl->GetYaxis()->SetTitle("N Clusters Fake");
\r
2760 nam = Form("%s_fake_cl",strN.Data());
\r
2761 tit = Form("%s Fake Clusters",strT.Data());
\r
2762 fHMCFakePatt = new TH1F(nam.Data(),tit.Data(),nlr,-0.5,nlr-0.5);
\r
2763 fHMCFakePatt->GetXaxis()->SetTitle("Fake clusters pattern");
\r
2764 fHMCFakePatt->Sumw2();
\r
2768 //____________________________________________________________________________
\r
2769 KMCTrackSummary::~KMCTrackSummary()
\r
2771 if (fHMCChi2) delete fHMCChi2;
\r
2772 if (fHMCSigDCARPhi) delete fHMCSigDCARPhi;
\r
2773 if (fHMCSigDCAZ) delete fHMCSigDCAZ;
\r
2774 if (fHMCSigPt) delete fHMCSigPt;
\r
2775 if (fHMCNCl) delete fHMCNCl;
\r
2776 if (fHMCFakePatt) delete fHMCFakePatt;
\r
2780 //____________________________________________________________________________
\r
2781 void KMCTrackSummary::Add(const KMCTrackSummary* src, double scl)
\r
2783 if (fHMCChi2) fHMCChi2->Add(src->fHMCChi2,scl);
\r
2784 if (fHMCSigDCARPhi) fHMCSigDCARPhi->Add(src->fHMCSigDCARPhi,scl);
\r
2785 if (fHMCSigDCAZ) fHMCSigDCAZ->Add(src->fHMCSigDCAZ,scl);
\r
2786 if (fHMCSigPt) fHMCSigPt->Add(src->fHMCSigPt,scl);
\r
2787 if (fHMCNCl) fHMCNCl->Add(src->fHMCNCl,scl);
\r
2788 if (fHMCFakePatt) fHMCFakePatt->Add(src->fHMCFakePatt,scl);
\r
2789 fCountAcc += src->fCountAcc;
\r
2790 fUpdCalls += src->fUpdCalls;
\r
2793 //____________________________________________________________________________
\r
2794 void KMCTrackSummary::PrependTNamed(TNamed* obj, const char* nm, const char* tit)
\r
2796 // prepend name, title by this prefix
\r
2797 TString name = obj->GetName(),title = obj->GetTitle();
\r
2798 name.Prepend(nm); title.Prepend(tit);
\r
2799 obj->SetNameTitle(name.Data(),title.Data());
\r
2803 //____________________________________________________________________________
\r
2804 void KMCTrackSummary::SetNamePrefix(const char* pref)
\r
2806 // prepend all names by this prefix
\r
2807 TString nmT,nmP = pref;
\r
2808 if (nmP.IsNull()) return;
\r
2809 nmP = Form("%s_",pref);
\r
2810 nmT = Form("%s ",pref);
\r
2811 PrependTNamed(this, nmP.Data(), nmT.Data());
\r
2812 PrependTNamed(fHMCChi2, nmP.Data(), nmT.Data());
\r
2813 PrependTNamed(fHMCSigDCARPhi, nmP.Data(), nmT.Data());
\r
2814 PrependTNamed(fHMCSigDCAZ, nmP.Data(), nmT.Data());
\r
2815 PrependTNamed(fHMCSigPt, nmP.Data(), nmT.Data());
\r
2816 PrependTNamed(fHMCNCl, nmP.Data(), nmT.Data());
\r
2817 PrependTNamed(fHMCFakePatt, nmP.Data(), nmT.Data());
\r
2821 //____________________________________________________________________________
\r
2822 Bool_t KMCTrackSummary::CheckTrack(KMCProbe* trc)
\r
2824 // check if the track satisfies to selection conditions
\r
2826 if (!trc) return kFALSE;
\r
2827 if (OutOfRange(trc->GetNITSHits(), fMinMaxClITS)) return kFALSE;
\r
2828 if (OutOfRange(trc->GetNFakeITSHits(), fMinMaxClITSFake)) return kFALSE;
\r
2829 if (OutOfRange(trc->GetNITSHits()-trc->GetNFakeITSHits(), fMinMaxClITSCorr)) return kFALSE;
\r
2831 // check layer patterns if requested
\r
2832 UInt_t patt = trc->GetHitsPatt();
\r
2833 UInt_t pattF = trc->GetFakesPatt();
\r
2834 UInt_t pattC = patt&(~pattF); // correct hits
\r
2835 // is there at least one hit in each of requested layer groups?
\r
2836 for (int ip=fPattITS.GetSize();ip--;) if (!CheckPattern(patt,fPattITS[ip])) return kFALSE;
\r
2837 // is there at least one fake in any of requested layer groups?
\r
2838 for (int ip=fPattITSFakeExcl.GetSize();ip--;) if (CheckPattern(pattF,fPattITSFakeExcl[ip])) return kFALSE;
\r
2839 // is there at least one correct hit in each of requested layer groups?
\r
2840 for (int ip=fPattITSCorr.GetSize();ip--;) if (!CheckPattern(pattC,fPattITSCorr[ip])) return kFALSE;
\r
2846 //____________________________________________________________________________
\r
2847 void KMCTrackSummary::AddTrack(KMCProbe* trc)
\r
2849 // fill track info
\r
2850 if (!CheckTrack(trc)) return;
\r
2852 fHMCChi2->Fill(trc->GetNormChi2(kFALSE));
\r
2853 fHMCSigDCARPhi->Fill(trc->GetY());
\r
2854 fHMCSigDCAZ->Fill(trc->GetZ());
\r
2855 if (fRefProbe.Pt()>0) fHMCSigPt->Fill( trc->Pt()/fRefProbe.Pt()-1.);
\r
2856 // printf("Pts: %.3f %.3f -> %.3f\n",trc->Pt(),fRefProbe.Pt(),trc->Pt()/fRefProbe.Pt()-1.);
\r
2857 // trc->Print("l");
\r
2858 // fRefProbe.Print("l");
\r
2859 fHMCNCl->Fill(trc->GetNITSHits(),trc->GetNFakeITSHits());
\r
2860 for (int i=trc->GetNITSLayers();i--;) if (trc->IsHitFake(i)) fHMCFakePatt->Fill(i);
\r
2864 //____________________________________________________________________________
\r
2865 UInt_t KMCTrackSummary::Bits(Bool_t l0, Bool_t l1, Bool_t l2, Bool_t l3, Bool_t l4, Bool_t l5, Bool_t l6, Bool_t l7,Bool_t l8, Bool_t l9,
\r
2866 Bool_t l10,Bool_t l11,Bool_t l12,Bool_t l13,Bool_t l14,Bool_t l15,Bool_t l16,Bool_t l17,Bool_t l18,Bool_t l19,
\r
2867 Bool_t l20,Bool_t l21,Bool_t l22,Bool_t l23,Bool_t l24,Bool_t l25,Bool_t l26,Bool_t l27,Bool_t l28,Bool_t l29,
\r
2868 Bool_t l30,Bool_t l31)
\r
2870 // create corresponding bit pattern
\r
2872 if (l0 ) patt |= (0x1<<0); if (l1 ) patt |= (0x1<<1); if (l2 ) patt |= (0x1<<2);
\r
2873 if (l3 ) patt |= (0x1<<3); if (l4 ) patt |= (0x1<<4); if (l5 ) patt |= (0x1<<5);
\r
2874 if (l6 ) patt |= (0x1<<6); if (l7 ) patt |= (0x1<<7); if (l8 ) patt |= (0x1<<8);
\r
2875 if (l9 ) patt |= (0x1<<9); if (l10) patt |= (0x1<<10); if (l11) patt |= (0x1<<11);
\r
2876 if (l12) patt |= (0x1<<12); if (l13) patt |= (0x1<<13); if (l14) patt |= (0x1<<14);
\r
2877 if (l15) patt |= (0x1<<15); if (l16) patt |= (0x1<<16); if (l17) patt |= (0x1<<17);
\r
2878 if (l18) patt |= (0x1<<18); if (l19) patt |= (0x1<<19); if (l20) patt |= (0x1<<20);
\r
2879 if (l21) patt |= (0x1<<21); if (l22) patt |= (0x1<<22); if (l23) patt |= (0x1<<23);
\r
2880 if (l24) patt |= (0x1<<24); if (l25) patt |= (0x1<<25); if (l26) patt |= (0x1<<26);
\r
2881 if (l27) patt |= (0x1<<27); if (l28) patt |= (0x1<<28); if (l29) patt |= (0x1<<29);
\r
2882 if (l30) patt |= (0x1<<30); if (l31) patt |= (0x1<<31);
\r
2886 //__________________________________________________________________________
\r
2887 void KMCTrackSummary::Print(Option_t* ) const
\r
2889 printf("%s: summary for track M=%5.3f pT: %6.3f eta: %.2f\n", GetName(),
\r
2890 fRefProbe.GetMass(),fRefProbe.Pt(), fRefProbe.Eta());
\r
2891 printf("Cuts on NCl ITS: Tot: %2d - %2d Fake: %2d - %2d Corr: %2d - %2d\n",
\r
2892 fMinMaxClITS[0],fMinMaxClITS[1],
\r
2893 fMinMaxClITSFake[0],fMinMaxClITSFake[1],
\r
2894 fMinMaxClITSCorr[0],fMinMaxClITSCorr[1]);
\r
2896 int nlr = fNITSLayers;
\r
2897 if (fPattITS.GetSize()) {
\r
2898 printf("Require at least 1 hit in groups: ");
\r
2899 printf("Hits obligatory in groups: ");
\r
2900 for (int i=fPattITS.GetSize();i--;) {
\r
2901 UInt_t pat = (UInt_t)fPattITS[i];
\r
2903 for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
\r
2909 if (fPattITSFakeExcl.GetSize()) {
\r
2910 printf("Fake hit forbidden in groups : ");
\r
2911 for (int i=fPattITSFakeExcl.GetSize();i--;) {
\r
2912 UInt_t pat = (UInt_t)fPattITSFakeExcl[i];
\r
2914 for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
\r
2920 if (fPattITSCorr.GetSize()) {
\r
2921 printf("Correct hit obligatory in groups: ");
\r
2922 for (int i=fPattITSCorr.GetSize();i--;) {
\r
2923 UInt_t pat = (UInt_t)fPattITSCorr[i];
\r
2925 for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
\r
2931 printf("Entries: Tot: %.4e Acc: %.4e -> Eff: %.4f+-%.4f %.2e KMC updates\n",fCountTot,fCountAcc,GetEff(),GetEffErr(),GetUpdCalls());
\r