1 #include "DetectorK.h"
\r
3 #include <TMatrixD.h>
\r
6 #include <TFormula.h>
\r
8 #include <TEllipse.h>
\r
10 #include <TGraphErrors.h>
\r
12 #include "AliExternalTrackParam.h"
\r
14 /***********************************************************
\r
16 Fast Simulation tool for Inner Tracker Systems
\r
18 original code of using the billoir technique was developed
\r
19 for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov
\r
20 http://rnc.lbl.gov/~jhthomas
\r
22 Changes by S. Rossegger -> see header file
\r
24 ***********************************************************/
\r
27 #define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
\r
29 #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
\r
30 #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm)
\r
31 #define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)
\r
32 // #define dNdEtaCent 2300//15000 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known)
\r
34 #define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
\r
35 #define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
\r
36 #define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
\r
37 #define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0)
\r
38 #define EfficiencySearchFlag 2 // Define search method:
\r
39 // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
\r
41 #define PionMass 0.139 // Mass of the Pion
\r
42 #define KaonMass 0.498 // Mass of the Kaon
\r
43 #define D0Mass 1.865 // Mass of the D0
\r
45 //TMatrixD *probKomb; // table for efficiency kombinatorics
\r
48 class CylLayerK : public TNamed {
\r
51 CylLayerK(char *name) : TNamed(name,name) {}
\r
53 Float_t GetRadius() const {return radius;}
\r
54 Float_t GetRadL() const {return radL;}
\r
55 Float_t GetPhiRes() const {return phiRes;}
\r
56 Float_t GetZRes() const {return zRes;}
\r
57 Float_t GetLayerEff() const {return eff;}
\r
59 // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
\r
60 Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes;
\r
64 ClassDef(CylLayerK,1);
\r
68 class ForwardLayer : public TNamed {
\r
70 ForwardLayer(char *name) : TNamed(name,name) {}
\r
72 Float_t GetZ() const {return zPos;}
\r
73 Float_t GetXRes() const {return xRes;}
\r
74 Float_t GetYRes() const {return yRes;}
\r
75 Float_t GetThickness() const {return thickness;}
\r
76 Float_t Getdensity() const {return density;}
\r
77 Float_t GetLayerEff() const {return eff;}
\r
79 // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
\r
80 Float_t zPos; Float_t xRes; Float_t yRes;
\r
87 ClassDef(ForwardLayer,1);
\r
92 DetectorK::DetectorK()
\r
93 : TNamed("test_detector","detector"),
\r
95 fNumberOfActiveLayers(0),
\r
96 fNumberOfActiveITSLayers(0),
\r
99 fIntegrationTime(0.02), // in ms
\r
100 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
\r
101 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
\r
102 fParticleMass(0.140), // Standard: pion mass
\r
103 fMaxRadiusSlowDet(10.),
\r
104 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
\r
105 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
\r
109 // default constructor
\r
111 // fLayers = new TObjArray();
\r
115 DetectorK::DetectorK(char *name, char *title)
\r
116 : TNamed(name,title),
\r
117 fNumberOfLayers(0),
\r
118 fNumberOfActiveLayers(0),
\r
119 fNumberOfActiveITSLayers(0),
\r
122 fIntegrationTime(0.02), // in ms
\r
123 fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
\r
124 fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
\r
125 fParticleMass(0.140), // Standard: pion mass
\r
126 fMaxRadiusSlowDet(10.),
\r
127 fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
\r
128 fAtLeastFake(1), // if at least x fakes, track is considered fake ...
\r
132 // default constructor, that set the name and title
\r
134 // fLayers = new TObjArray();
\r
136 DetectorK::~DetectorK() { //
\r
137 // virtual destructor
\r
142 void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) {
\r
144 // Add additional layer to the list of layers (ordered by radius)
\r
147 CylLayerK *newLayer = (CylLayerK*) fLayers.FindObject(name);
\r
150 newLayer = new CylLayerK(name);
\r
151 newLayer->radius = radius;
\r
152 newLayer->radL = radL;
\r
153 newLayer->phiRes = phiRes;
\r
154 newLayer->zRes = zRes;
\r
155 newLayer->eff = eff;
\r
157 if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS)
\r
158 newLayer->isDead = kTRUE;
\r
160 newLayer->isDead = kFALSE;
\r
162 if (fLayers.GetEntries()==0)
\r
163 fLayers.Add(newLayer);
\r
166 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
167 CylLayerK *l = (CylLayerK*)fLayers.At(i);
\r
168 if (radius<l->radius) {
\r
169 fLayers.AddBefore(l,newLayer);
\r
172 if (radius>l->radius && (i+1)==fLayers.GetEntries() ) {
\r
173 // even bigger then last one
\r
174 fLayers.Add(newLayer);
\r
179 fNumberOfLayers += 1;
\r
180 if (!(newLayer->isDead)) {
\r
181 fNumberOfActiveLayers += 1;
\r
182 TString lname(newLayer->GetName());
\r
183 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
\r
188 printf("Layer with the name %s does already exist\n",name);
\r
194 void DetectorK::KillLayer(char *name) {
\r
196 // Marks layer as dead. Contribution only by Material Budget
\r
199 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
201 printf("Layer %s not found - cannot mark as dead\n",name);
\r
203 tmp->phiRes = 999999;
\r
204 tmp->zRes = 999999;
\r
205 if (!(tmp->isDead)) {
\r
206 tmp->isDead = kTRUE;
\r
207 fNumberOfActiveLayers -= 1;
\r
208 TString lname(tmp->GetName());
\r
209 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
\r
214 void DetectorK::SetRadius(char *name, Float_t radius) {
\r
216 // Set layer radius [cm]
\r
219 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
223 printf("Layer %s not found - cannot set radius\n",name);
\r
226 Float_t tmpRadL = tmp->radL;
\r
227 Float_t tmpPhiRes = tmp->phiRes;
\r
228 Float_t tmpZRes = tmp->zRes;
\r
230 RemoveLayer(name); // so that the ordering is correct
\r
231 AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes);
\r
235 Float_t DetectorK::GetRadius(char *name) {
\r
237 // Return layer radius [cm]
\r
240 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
242 printf("Layer %s not found - cannot get radius\n",name);
\r
244 return tmp->radius;
\r
249 void DetectorK::SetRadiationLength(char *name, Float_t radL) {
\r
251 // Set layer material [cm]
\r
254 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
256 printf("Layer %s not found - cannot set layer material\n",name);
\r
262 Float_t DetectorK::GetRadiationLength(char *name) {
\r
264 // Return layer radius [cm]
\r
267 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
269 printf("Layer %s not found - cannot get layer material\n",name);
\r
277 void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
\r
279 // Set layer resolution in [cm]
\r
282 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
284 printf("Layer %s not found - cannot set resolution\n",name);
\r
287 Bool_t wasDead = tmp->isDead;
\r
289 tmp->phiRes = phiRes;
\r
291 TString lname(tmp->GetName());
\r
293 if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {
\r
294 tmp->isDead = kTRUE;
\r
296 fNumberOfActiveLayers -= 1;
\r
297 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
\r
300 tmp->isDead = kFALSE;
\r
302 fNumberOfActiveLayers += 1;
\r
303 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
\r
311 Float_t DetectorK::GetResolution(char *name, Int_t axis) {
\r
313 // Return layer resolution in [cm]
\r
314 // axis = 0: resolution in rphi
\r
315 // axis = 1: resolution in z
\r
318 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
320 printf("Layer %s not found - cannot get resolution\n",name);
\r
322 if (axis==0) return tmp->phiRes;
\r
323 if (axis==1) return tmp->zRes;
\r
324 printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
\r
329 void DetectorK::SetLayerEfficiency(char *name, Float_t eff) {
\r
331 // Set layer efficnecy (prop that his is missed within this layer)
\r
334 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
336 printf("Layer %s not found - cannot set layer efficiency\n",name);
\r
342 Float_t DetectorK::GetLayerEfficiency(char *name) {
\r
344 // Get layer efficnecy (prop that his is missed within this layer)
\r
347 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
349 printf("Layer %s not found - cannot get layer efficneicy\n",name);
\r
357 void DetectorK::RemoveLayer(char *name) {
\r
359 // Removes a layer from the list
\r
362 CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
\r
364 printf("Layer %s not found - cannot remove it\n",name);
\r
366 Bool_t wasDead = tmp->isDead;
\r
367 fLayers.Remove(tmp);
\r
368 fNumberOfLayers -= 1;
\r
370 fNumberOfActiveLayers -= 1;
\r
371 TString lname(tmp->GetName());
\r
372 if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
\r
379 void DetectorK::PrintLayout() {
\r
381 // Prints the detector layout
\r
384 printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
\r
386 if (fLayers.GetEntries()>0)
\r
387 printf(" Name \t\t r [cm] \t X0 \t phi & z res [um]\n");
\r
389 CylLayerK *tmp = 0;
\r
390 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
391 tmp = (CylLayerK*)fLayers.At(i);
\r
393 // don't print all the tpc layers
\r
394 TString name(tmp->GetName());
\r
395 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
\r
397 printf("%d. %s \t %03.2f \t%1.4f\t ",i,
\r
398 tmp->GetName(), tmp->radius, tmp->radL);
\r
399 if (tmp->phiRes==RIDICULOUS)
\r
402 printf("%3.0f ",tmp->phiRes*10000);
\r
403 if (tmp->zRes==RIDICULOUS)
\r
406 printf("%3.0f\n",tmp->zRes*10000);
\r
410 void DetectorK::PlotLayout(Int_t plotDead) {
\r
412 // Plots the detector layout in Front view
\r
415 Double_t x0=0, y0=0;
\r
417 TGraphErrors *gr = new TGraphErrors();
\r
418 gr->SetPoint(0,0,0);
\r
419 CylLayerK *lastLayer = (CylLayerK*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->radius;
\r
420 gr->SetPointError(0,maxRad,maxRad);
\r
424 CylLayerK *tmp = 0;
\r
425 for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
\r
426 tmp = (CylLayerK*)fLayers.At(i);
\r
429 Double_t txtpos = tmp->radius;
\r
430 if ((tmp->isDead)) txtpos*=-1; //
\r
431 TText *txt = new TText(x0,txtpos,tmp->GetName());
\r
432 txt->SetTextSizePixels(5); txt->SetTextAlign(21);
\r
433 if (!tmp->isDead || plotDead) txt->Draw();
\r
435 TEllipse *layEl = new TEllipse(x0,y0,tmp->radius);
\r
436 // layEl->SetFillColor(5);
\r
437 layEl->SetFillStyle(5001);
\r
438 layEl->SetLineStyle(tmp->isDead+1); // dashed if not active
\r
439 layEl->SetLineColor(4);
\r
440 TString name(tmp->GetName());
\r
441 if (!tmp->isDead) layEl->SetLineWidth(2);
\r
442 if (name.Contains("tpc") ) layEl->SetLineColor(29);
\r
444 if (!tmp->isDead || plotDead) layEl->Draw();
\r
452 void DetectorK::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
\r
454 // Emulates the TPC
\r
456 // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
\r
459 AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage
\r
461 // % Radiation Lengths ... Average per TPC row (i.e. total/159 )
\r
462 Float_t radLPerRow = 0.000036;
\r
464 Float_t tpcInnerRadialPitch = 0.75 ; // cm
\r
465 Float_t tpcMiddleRadialPitch = 1.0 ; // cm
\r
466 Float_t tpcOuterRadialPitch = 1.5 ; // cm
\r
467 // Float_t tpcInnerPadWidth = 0.4 ; // cm
\r
468 // Float_t tpcMiddlePadWidth = 0.6 ; // cm
\r
469 // Float_t tpcOuterPadWidth = 0.6 ; // cm
\r
470 Float_t innerRows = 63 ;
\r
471 Float_t middleRows = 64 ;
\r
472 Float_t outerRows = 32 ;
\r
473 Float_t tpcRows = (innerRows + middleRows + outerRows) ;
\r
474 Float_t rowOneRadius = 85.2 ; // cm
\r
475 Float_t row64Radius = 135.1 ; // cm
\r
476 Float_t row128Radius = 199.2 ; // cm
\r
478 for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
\r
480 Float_t rowRadius =0;
\r
482 rowRadius = rowOneRadius + k*tpcInnerRadialPitch ;
\r
483 else if ( k>=innerRows && k<(innerRows+middleRows) )
\r
484 rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
\r
485 else if (k>=(innerRows+middleRows) && k<tpcRows )
\r
486 rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
\r
489 AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean);
\r
491 AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row
\r
498 void DetectorK::RemoveTPC() {
\r
500 // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
\r
501 CylLayerK *tmp = 0;
\r
502 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
503 tmp = (CylLayerK*)fLayers.At(i);
\r
504 TString name(tmp->GetName());
\r
505 if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
\r
507 RemoveLayer((char*)"IFC");
\r
512 Double_t DetectorK::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const
\r
515 // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
\r
518 Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ;
\r
519 Double_t theta = 0.0 ; // Momentum and mass in GeV
\r
520 // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
\r
521 if ( radLength > 0 ) theta = 0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ;
\r
526 Double_t DetectorK::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
528 // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
\r
529 // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
\r
530 // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
\r
531 Double_t sx, sy, goodHit ;
\r
532 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
\r
533 sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ;
\r
534 goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ;
\r
535 return ( goodHit ) ;
\r
539 Double_t DetectorK::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
541 // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
\r
542 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
\r
543 Double_t sx, goodHit ;
\r
544 sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
\r
545 goodHit = 1./(1+sx) ;
\r
546 return ( goodHit ) ;
\r
549 Double_t DetectorK::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
551 // Based on work by Ruben Shahoyen
\r
552 // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
\r
553 // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
\r
554 // Following is correct for 2 DOF
\r
556 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
\r
557 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
\r
558 Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));
\r
559 return ( goodHit ) ;
\r
562 Double_t DetectorK::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
\r
564 // Based on work by Ruben Shahoyen
\r
565 // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)
\r
567 Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
\r
568 Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
\r
569 Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));
\r
570 return ( nullHit ) ;
\r
573 Double_t DetectorK::HitDensity ( Double_t radius )
\r
575 // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
\r
576 // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
\r
577 // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
\r
578 // Note that this function assumes we are working in CM and CM**2 [not meters].
\r
579 // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
\r
581 // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors
\r
582 // and only fast detectors reside outside this radius.
\r
583 Double_t arealDensity = 0 ;
\r
585 if ( radius > fMaxRadiusSlowDet )
\r
587 arealDensity = OneEventHitDensity(fdNdEtaCent,radius) ; // Fast detectors see central collision density (only)
\r
588 arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background
\r
591 if (radius < fMaxRadiusSlowDet )
\r
592 { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
\r
593 arealDensity = OneEventHitDensity(fdNdEtaCent,radius)
\r
594 + IntegratedHitDensity(dNdEtaMinB,radius)
\r
595 + UpcHitDensity(radius) ;
\r
596 arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
\r
597 // Increase density due to background
\r
600 return ( arealDensity ) ;
\r
604 double DetectorK::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
\r
606 // This is for one event at the vertex. No smearing.
\r
608 double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
\r
609 double tg = TMath::Tan(2*TMath::ATan(TMath::Exp(-fAvgRapidity)));
\r
610 den = den/TMath::Sqrt(1 + 1/(tg*tg));
\r
612 // double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
\r
613 // note: surface of sphere is '4*pi*r^2'
\r
614 // surface of cylinder is '2*pi*r* h'
\r
622 double DetectorK::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
\r
624 // The integral of minBias events smeared over a gaussian vertex distribution.
\r
625 // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
\r
627 Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
\r
628 Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
\r
630 // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
\r
631 if ( den < OneEventHitDensity(multiplicity,radius) ) den = OneEventHitDensity(multiplicity,radius) ;
\r
637 double DetectorK::UpcHitDensity(Double_t radius)
\r
639 // QED electrons ...
\r
641 Double_t mUPCelectrons ; ;
\r
642 // mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
\r
643 mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB; // Fit to 'Rossegger,Sadovsky'-Alice simulation
\r
644 if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R
\r
645 mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate
\r
646 mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC
\r
648 return mUPCelectrons ;
\r
652 double DetectorK::Dist(double z, double r)
\r
654 // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution
\r
655 // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
\r
656 // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
\r
657 Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration
\r
658 Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work
\r
659 double dist = 0.0 ;
\r
660 double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm
\r
661 double z0 = 0.0 ; //cm
\r
662 for(int i=0; i<nsteps; i++){
\r
663 if ( i == nsteps-1 ) index = 1 ;
\r
664 z0 = -4*SigmaD + i*dz0 ;
\r
665 dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
\r
666 (1/sqrt((z-z0)*(z-z0) + r*r)) ;
\r
667 if ( index != 4 ) index = 4; else index = 2 ;
\r
672 #define PZero 0.861 // Momentum of back to back decay particles in the CM frame
\r
673 #define EPiZero 0.872 // Energy of the pion from a D0 decay at rest
\r
674 #define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest
\r
676 Double_t DetectorK::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const {
\r
677 // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi.
\r
679 Double_t const1 = pt / D0Mass ;
\r
680 Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
\r
681 Double_t sum, ptPi, ptK ;
\r
682 Double_t effp, effk ;
\r
685 for ( Int_t k = 0 ; k < 360 ; k++ ) {
\r
687 Double_t theta = k * TMath::Pi() / 180. ;
\r
689 ptPi = TMath::Sqrt(
\r
690 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
\r
691 const1*const1*EPiZero*EPiZero -
\r
692 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
\r
693 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
\r
696 ptK = TMath::Sqrt(
\r
697 PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
\r
698 const1*const1*EKZero*EKZero +
\r
699 2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
\r
700 PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
\r
703 // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
\r
704 Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
\r
705 Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
\r
707 if ( pionindex >= kNptBins ) pionindex = 399 ;
\r
708 if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ;
\r
709 if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
\r
710 if ( effp < 0 ) effp = 0 ;
\r
712 if ( kaonindex >= kNptBins ) kaonindex = 399 ;
\r
713 if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ;
\r
714 if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
\r
715 if ( effk < 0 ) effk = 0 ;
\r
717 // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
\r
719 sum += effp * effk ;
\r
723 Double_t mean =sum/360;
\r
729 void DetectorK::SolveDOFminusOneAverage() {
\r
731 // Short study to address "# layers-1 efficiencies"
\r
732 // saves the means in the according arrays
\r
733 // Note: Obviously, does not work for the Telescope equation
\r
736 Double_t fMomentumResM[kNptBins], fResolutionRPhiM[kNptBins], fResolutionZM[kNptBins];
\r
737 Double_t efficiencyM[3][kNptBins];
\r
738 for (Int_t i=0; i<kNptBins; i++) {
\r
739 fMomentumResM[i] = 0; // Momentum resolution
\r
740 fResolutionRPhiM[i] = 0; // Resolution in R
\r
741 fResolutionZM[i] = 0; // Resolution in Z
\r
742 for (Int_t part=0; part<3; part++)
\r
743 efficiencyM[part][i] = 0; // efficiencies
\r
746 // loop over active layers in ITS (remove 1 by 1)
\r
747 Int_t nITSLayers = 0;
\r
748 CylLayerK *layer =0;
\r
749 for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) {
\r
750 layer = (CylLayerK*)fLayers.At(j);
\r
751 TString name(layer->GetName());
\r
752 if (name.Contains("tpc")) continue;
\r
753 if (!(layer->isDead)) {
\r
756 printf("Kill Layer %s\n",name.Data());
\r
757 Double_t rRes = GetResolution((char*)name.Data(),0);
\r
758 Double_t zRes = GetResolution((char*)name.Data(),1);
\r
759 KillLayer((char*)name.Data());
\r
761 SolveViaBilloir(1,0);
\r
763 // produce sum for the mean calculation
\r
764 for (Int_t i=0; i<kNptBins; i++) {
\r
765 fMomentumResM[i] += fMomentumRes[i]; // Momentum resolution
\r
766 fResolutionRPhiM[i] += fResolutionRPhi[i]; // Resolution in R
\r
767 fResolutionZM[i] += fResolutionZ[i]; // Resolution in Z
\r
768 for (Int_t part=0; part<3; part++)
\r
769 efficiencyM[part][i] += fEfficiency[part][i]; // efficiencies
\r
772 // "Restore" layer ...
\r
773 SetResolution((char*)name.Data(),rRes,zRes);
\r
778 // save means in "std. Arrays"
\r
779 for (Int_t i=0; i<kNptBins; i++) {
\r
780 fMomentumRes[i] = fMomentumResM[i]/nITSLayers; // Momentum resolution
\r
781 fResolutionRPhi[i] = fResolutionRPhiM[i]/nITSLayers; // Resolution in R
\r
782 fResolutionZ[i] = fResolutionZM[i]/nITSLayers; // Resolution in Z
\r
783 for (Int_t part=0; part<3; part++)
\r
784 fEfficiency[part][i] = efficiencyM[part][i]/nITSLayers; // efficiencies
\r
790 void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt) {
\r
792 // Solves the current geometry with the Billoir technique
\r
793 // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. )
\r
794 // ABOVE IS OBSOLETE -> NOW, its uses the Aliroot Kalman technique
\r
797 static AliExternalTrackParam probTr; // track to propagate
\r
798 probTr.SetUseLogTermMS(kTRUE);
\r
801 Int_t nPt = kNptBins;
\r
803 for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {
\r
804 for (Int_t j=0; j<nPt; j++) {
\r
805 fDetPointRes[i][j] = RIDICULOUS;
\r
806 fDetPointZRes[i][j] = RIDICULOUS;
\r
807 fTransMomenta[i] =0;
\r
808 fMomentumRes[i] =0;
\r
809 fResolutionRPhi[i] =0;
\r
813 if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt'
\r
818 // Calculate track parameters using Billoirs method of matrices
\r
820 Double_t pt,tgl, pz, lambda, deltaPoverP ;
\r
823 Int_t printOnce = 1 ;
\r
825 mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k
\r
827 mass[2] = fParticleMass; // third loop
\r
830 if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode
\r
835 // Prepare Probability Kombinations
\r
836 Int_t nLayer = fNumberOfActiveITSLayers;
\r
837 Int_t base = 3; // null, fake, correct
\r
839 Int_t komb = TMath::Power(base,nLayer);
\r
841 TMatrixD probLay(base,fNumberOfActiveITSLayers);
\r
842 TMatrixD probKomb(komb,nLayer);
\r
843 for (Int_t num=0; num<komb; num++) {
\r
844 for (Int_t l=nLayer-1; l>=0; l--) {
\r
845 Int_t pow = ((Int_t)TMath::Power(base,l+1));
\r
846 probKomb(num,nLayer-1-l)=(num%pow)/((Int_t)TMath::Power(base,l));
\r
850 for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) {
\r
852 // PseudoRapidity OK, used as an angle
\r
853 lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ;
\r
856 for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop
\r
858 CylLayerK *last = (CylLayerK*) fLayers.At((fLayers.GetEntries()-1));
\r
860 // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV
\r
861 Double_t bigRad = last->radius/2 ; // min. pt which the algorithm below could handle
\r
862 //if (bigRad<61) bigRad=61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) )
\r
863 fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ;
\r
864 if (!allPt) { // just 3 points around meanPt
\r
865 fTransMomenta[i] = meanPt-0.01+(Double_t)(i)*0.01;
\r
868 // New from here ................
\r
870 // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
\r
871 // These are the EndPoint values for y, z, a, b, and d
\r
872 double bGauss = fBField*10; // field in kgauss
\r
873 pt = fTransMomenta[i]; // GeV/c
\r
874 tgl = TMath::Tan(lambda); // dip
\r
875 charge = -1; // Assume an electron
\r
876 pz = pt * TMath::Tan(lambda) ; // GeV/
\r
877 enum {kY,kZ,kSnp,kTgl,kPtI}; // track parameter aliases
\r
878 enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
\r
881 double *trPars = (double*)probTr.GetParameter();
\r
882 double *trCov = (double*)probTr.GetCovariance();
\r
883 trPars[kY] = 0; // start from Y = 0
\r
884 trPars[kZ] = 0; // Z = 0
\r
885 trPars[kSnp] = 0; // track along X axis at the vertex
\r
886 trPars[kTgl] = TMath::Tan(lambda); // dip
\r
887 trPars[kPtI] = charge/pt; // q/pt
\r
889 // put tiny errors to propagate to the outer radius
\r
890 trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;
\r
892 if (!GetXatLabR(&probTr, last->radius ,xR,bGauss,1)) {
\r
893 printf("Track with pt=%f cannot reach radius %f\n",pt,last->radius);
\r
896 probTr.PropagateTo(xR, bGauss); // bring track to outer layer
\r
897 // reset cov.matrix
\r
898 const double kLargeErr2Coord = 50*50;
\r
899 const double kLargeErr2Dir = 0.6*0.6;
\r
900 const double kLargeErr2PtI = 0.5*0.5;
\r
901 for (int ic=15;ic--;) trCov[ic] = 0.;
\r
902 trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
\r
903 trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
\r
904 trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
\r
906 // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);
\r
908 // Set Detector-Efficiency Storage area to unity
\r
909 fEfficiency[massloop][i] = 1.0 ;
\r
911 // Back-propagate the covariance matrix along the track.
\r
913 CylLayerK *layer = 0;
\r
915 // find last "active layer" - start tracking at the last active layer
\r
916 Int_t lastActiveLayer = 0;
\r
917 for (Int_t j=fLayers.GetEntries(); j--;) {
\r
918 layer = (CylLayerK*)fLayers.At(j);
\r
919 if (!(layer->isDead)) { // is alive
\r
920 lastActiveLayer = j;
\r
925 for (Int_t j=lastActiveLayer; j--;) { // Layer loop
\r
927 layer = (CylLayerK*)fLayers.At(j);
\r
928 TString name(layer->GetName());
\r
929 Bool_t isVertex = name.Contains("vertex");
\r
931 if (!GetXatLabR(&probTr, layer->radius ,xR,bGauss)) {
\r
932 printf("Track with pt=%f cannot reach radius %f. This should not happen here\n",pt,layer->radius);
\r
936 probTr.PropagateTo(xR, bGauss); // propagate to this layer
\r
938 // rotate to frame with X axis normal to the surface
\r
941 probTr.GetXYZ(pos); // lab position
\r
942 double phi = TMath::ATan2(pos[1],pos[0]);
\r
943 if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
\r
944 if (!probTr.Rotate(phi)) {
\r
945 printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f\n",
\r
946 phi,layer->radius,pos[0],pos[1],pos[2]);
\r
951 // save resolutions at this layer
\r
952 fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
\r
953 fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
\r
954 //printf(">> L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
\r
957 if (isVertex) continue;
\r
959 // create fake measurement with the errors assigned to the layer
\r
960 // account for the measurement there
\r
961 double meas[2] = {probTr.GetY(),probTr.GetZ()};
\r
962 double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
\r
964 if (!probTr.Update(meas,measErr2)) {
\r
965 printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
\r
966 meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
\r
971 // correct for materials of this layer
\r
972 // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
\r
973 if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
\r
974 printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
\r
981 // Pattern recognition is done .... save values like vertex resolution etc.
\r
983 // Convert the Convariance matrix parameters into physical quantities
\r
984 // The results are propogated to the previous point but *do not* include the measurement at that point.
\r
985 // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P(); // Absolute magnitude so ignore charge
\r
986 deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P();
\r
987 fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent
\r
988 fResolutionRPhi[i] = TMath::Sqrt( probTr.GetSigmaY2() ) * 1.e4; // result in microns
\r
989 fResolutionZ[i] = TMath::Sqrt( probTr.GetSigmaZ2() ) * 1.e4; // result in microns
\r
990 // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
\r
993 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
994 printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;
\r
995 if (fAtLeastCorr != -1) printf("Number of combinatorics for probabilities: %d\n",komb);
\r
996 printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);
\r
997 printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ;
\r
1001 // print out and efficiency calculation
\r
1002 Int_t iLayActive=0;
\r
1003 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
\r
1005 layer = (CylLayerK*)fLayers.At(j);
\r
1007 // Convert to Meters, Tesla, and GeV
\r
1008 Float_t radius = layer->radius /100;
\r
1009 Float_t phiRes = layer->phiRes /100;
\r
1010 Float_t zRes = layer->zRes /100;
\r
1011 Float_t radLength = layer->radL;
\r
1012 Float_t leff = layer->eff; // basic layer efficiency
\r
1013 Bool_t isDead = layer->isDead;
\r
1016 if ( (!isDead && radLength >0) ) {
\r
1018 Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
\r
1019 phiRes * phiRes ) * 100. ; // work in cm
\r
1020 Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
\r
1021 zRes * zRes ) * 100. ; // work in cm
\r
1024 Double_t layerEfficiency = 0;
\r
1025 if ( EfficiencySearchFlag == 0 )
\r
1026 layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
\r
1027 else if ( EfficiencySearchFlag == 1 )
\r
1028 layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
\r
1029 else if ( EfficiencySearchFlag == 2 )
\r
1030 layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
\r
1032 TString name(layer->GetName());
\r
1033 if (!name.Contains("tpc")) {
\r
1034 probLay(2,iLayActive)= layerEfficiency ; // Pcorr
\r
1035 probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
\r
1036 probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
\r
1039 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
\r
1041 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
1044 printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
\r
1045 layer->GetName(), radius*100, radLength,
\r
1046 fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
\r
1047 phiRes*1.e6, zRes*1.e6,
\r
1048 HitDensity(radius*100)) ;
\r
1049 if (!name.Contains("tpc"))
\r
1050 printf("%10.3f\n", layerEfficiency);
\r
1055 if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;
\r
1060 if (fAtLeastCorr != -1) {
\r
1061 // Calculate probabilities from Kombinatorics tree ...
\r
1062 Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
\r
1063 fEfficiency[massloop][i] = probs[0]; // efficiency
\r
1064 fFake[massloop][i] = probs[1]; // fake
\r
1069 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {
\r
1070 printf("%s:\t ----- ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6);
\r
1074 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
1075 if (fNumberOfActiveLayers >=15) printOnce = 0 ;
\r
1082 if (fNumberOfActiveLayers <15 ) {
\r
1086 // BACKWORD TRACKING +++++++++++++++++
\r
1087 // number of layers is quite low ... efficiency calculation was probably nonsense
\r
1088 // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates"
\r
1090 // For below, see paper, NIM A262 (1987) p.444, eqs.12.
\r
1091 // Equivalently, one can simply combine the forward and backward estimates. Assuming
\r
1092 // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can
\r
1093 // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1).
\r
1094 // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.
\r
1096 Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot
\r
1098 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
1099 printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n",
\r
1100 fNumberOfActiveLayers);
\r
1103 // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam)
\r
1104 // mIstar.UnitMatrix(); // start with unity
\r
1105 if (doLikeAliRoot) {
\r
1106 probTr.ResetCovariance(10);
\r
1108 // cannot do complete reset, set to very large errors
\r
1109 trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
\r
1110 trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
\r
1111 trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
\r
1112 // cout<<pt<<": "<<kLargeErr2Coord<<" "<<kLargeErr2Dir<<" "<<kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]<<endl;
\r
1114 // Clean up and storing of "forward estimates"
\r
1115 Double_t detPointResForw[kMaxNumberOfDetectors][kNptBins], detPointZResForw[kMaxNumberOfDetectors][kNptBins] ;
\r
1116 for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {
\r
1117 for (Int_t l=0; l<nPt; l++) {
\r
1118 detPointResForw[k][l] = fDetPointRes[k][l];
\r
1119 if (!doLikeAliRoot) fDetPointRes[k][l] = RIDICULOUS;
\r
1120 detPointZResForw[k][l] = fDetPointZRes[k][l];
\r
1121 if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS;
\r
1125 // find first "active layer" - start tracking at the first active layer
\r
1126 Int_t firstActiveLayer = 0;
\r
1127 for (Int_t j=0; j<fLayers.GetEntries(); j++) {
\r
1128 layer = (CylLayerK*)fLayers.At(j);
\r
1129 if (!(layer->isDead)) { // is alive
\r
1130 firstActiveLayer = j;
\r
1135 for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()); j++) { // Layer loop
\r
1137 layer = (CylLayerK*)fLayers.At(j);
\r
1138 // CylLayerK *nextlayer = (CylLayerK*)fLayers.At(j+1);
\r
1140 TString name(layer->GetName());
\r
1141 Bool_t isVertex = name.Contains("vertex");
\r
1143 if (!GetXatLabR(&probTr, layer->radius ,xR,bGauss)) {
\r
1144 printf("Track with pt=%f cannot reach radius %f. This should not happen here\n",pt,layer->radius);
\r
1148 probTr.PropagateTo(xR, bGauss); // propagate to this layer
\r
1150 // rotate to frame with X axis normal to the surface
\r
1152 probTr.GetXYZ(pos); // lab position
\r
1153 double phi = TMath::ATan2(pos[1],pos[0]);
\r
1154 if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
\r
1155 if (!probTr.Rotate(phi)) {
\r
1156 printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f\n",
\r
1157 phi,layer->radius,pos[0],pos[1],pos[2]);
\r
1163 fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
\r
1164 fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
\r
1166 //printf("<< L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
\r
1167 // create fake measurement with the errors assigned to the layer
\r
1168 // account for the measurement there
\r
1169 if (isVertex) continue;
\r
1170 double meas[2] = {probTr.GetY(),probTr.GetZ()};
\r
1171 double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
\r
1173 if (!probTr.Update(meas,measErr2)) {
\r
1174 printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
\r
1175 meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
\r
1179 // correct for materials of this layer
\r
1180 // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
\r
1181 if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
\r
1182 printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
\r
1188 // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!!
\r
1190 // also update the values for the track position ??????
\r
1192 // Pattern recognition is done .... save values like vertex resolution etc.
\r
1194 // Invert the Matrix to recover the convariance matrix
\r
1196 // Convert the Convariance matrix parameters into physical quantities
\r
1197 // The results are propogated to the previous point but *do not* include the measurement at that point.
\r
1198 deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge
\r
1199 fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent
\r
1200 fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns
\r
1201 fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns
\r
1202 // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
\r
1205 // Weighted combination of the forward and backward estimates
\r
1206 if (!doLikeAliRoot) {
\r
1207 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {
\r
1208 fDetPointRes[j][i] = 1/(1/detPointResForw[j][i] + 1/fDetPointRes[j][i]);
\r
1209 fDetPointZRes[j][i] = 1/(1/detPointZResForw[j][i] + 1/fDetPointZRes[j][i]);
\r
1212 // Set Detector-Efficiency Storage area to unity
\r
1213 fEfficiency[massloop][i] = 1.0 ;
\r
1215 // print out and efficiency calculation
\r
1217 for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
\r
1219 layer = (CylLayerK*)fLayers.At(j);
\r
1221 // Convert to Meters, Tesla, and GeV
\r
1222 Float_t radius = layer->radius /100;
\r
1223 Float_t phiRes = layer->phiRes /100;
\r
1224 Float_t zRes = layer->zRes /100;
\r
1225 Float_t radLength = layer->radL;
\r
1226 Float_t leff = layer->eff;
\r
1227 Bool_t isDead = layer->isDead;
\r
1230 if ( (!isDead && radLength >0) ) {
\r
1231 Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
\r
1232 phiRes * phiRes ) * 100. ; // work in cm
\r
1233 Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
\r
1234 zRes * zRes ) * 100. ; // work in cm
\r
1236 Double_t layerEfficiency = 0;
\r
1237 if ( EfficiencySearchFlag == 0 )
\r
1238 layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
\r
1239 else if ( EfficiencySearchFlag == 1 )
\r
1240 layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
\r
1241 else if ( EfficiencySearchFlag == 2 )
\r
1242 layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
\r
1244 TString name(layer->GetName());
\r
1245 if (!name.Contains("tpc")) {
\r
1246 probLay(2,iLayActive)= layerEfficiency ; // Pcorr
\r
1247 probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
\r
1248 probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
\r
1251 if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
\r
1253 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
1256 printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
\r
1257 layer->GetName(), radius*100, radLength,
\r
1258 fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
\r
1259 phiRes*1.e6, zRes*1.e6,
\r
1260 HitDensity(radius*100)) ;
\r
1261 if (!name.Contains("tpc"))
\r
1262 printf("%10.3f\n", layerEfficiency);
\r
1267 if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;
\r
1270 if (fAtLeastCorr != -1) {
\r
1271 // Calculate probabilities from Kombinatorics tree ...
\r
1272 Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
\r
1273 fEfficiency[massloop][i] = probs[0]; // efficiency
\r
1274 fFake[massloop][i] = probs[1]; // fake
\r
1277 if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
\r
1285 probTr.SetUseLogTermMS(kFALSE); // Reset of MS term usage to avoid problems since its static
\r
1291 TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
\r
1293 // returns the momentum resolution
\r
1296 TGraph *graph = new TGraph(kNptBins, fTransMomenta, fMomentumRes);
\r
1297 graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
\r
1298 // graph->GetXaxis()->SetRangeUser(0.,5.0) ;
\r
1299 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1300 graph->GetXaxis()->CenterTitle();
\r
1301 graph->GetXaxis()->SetNoExponent(1) ;
\r
1302 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1303 graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;
\r
1304 graph->GetYaxis()->CenterTitle();
\r
1306 graph->SetMaximum(20) ;
\r
1307 graph->SetMinimum(0.1) ;
\r
1308 graph->SetLineColor(color);
\r
1309 graph->SetMarkerColor(color);
\r
1310 graph->SetLineWidth(linewidth);
\r
1316 TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
\r
1318 // Returns the pointing resolution
\r
1319 // axis = 0 ... rphi pointing resolution
\r
1320 // axis = 1 ... z pointing resolution
\r
1323 TGraph * graph = 0;
\r
1326 graph = new TGraph ( kNptBins, fTransMomenta, fResolutionRPhi ) ;
\r
1327 graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
\r
1328 graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
\r
1330 graph = new TGraph ( kNptBins, fTransMomenta, fResolutionZ ) ;
\r
1331 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
\r
1332 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
\r
1335 graph->SetMinimum(1) ;
\r
1336 graph->SetMaximum(300.1) ;
\r
1337 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1338 graph->GetXaxis()->CenterTitle();
\r
1339 graph->GetXaxis()->SetNoExponent(1) ;
\r
1340 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1341 graph->GetYaxis()->CenterTitle();
\r
1343 graph->SetLineWidth(linewidth);
\r
1344 graph->SetLineColor(color);
\r
1345 graph->SetMarkerColor(color);
\r
1352 TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
\r
1354 // returns the Pointing resolution (accoring to Telescope equation)
\r
1355 // axis =0 ... in rphi
\r
1356 // axis =1 ... in z
\r
1359 Double_t resolution[kNptBins];
\r
1361 Double_t layerResolution[2];
\r
1362 Double_t layerRadius[2];
\r
1363 Double_t layerThickness[2];
\r
1365 Int_t count =0; // search two first active layers
\r
1366 printf("Telescope equation for layers: ");
\r
1367 for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
\r
1368 CylLayerK *l = (CylLayerK*)fLayers.At(i);
\r
1369 if (!l->isDead && l->radius>0) {
\r
1370 layerRadius[count] = l->radius;
\r
1371 layerThickness[count] = l->radL;
\r
1373 layerResolution[count] = l->phiRes;
\r
1375 layerResolution[count] = l->zRes;
\r
1377 printf("%s, ",l->GetName());
\r
1380 if (count>=2) break;
\r
1384 Double_t pt, momentum, thickness,aMCS ;
\r
1385 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1387 for ( Int_t i = 0 ; i < kNptBins ; i++ ) {
\r
1388 // Reference data as if first two layers were acting all alone
\r
1389 pt = fTransMomenta[i] ;
\r
1390 momentum = pt / TMath::Cos(lambda) ; // Total momentum
\r
1391 resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
\r
1392 + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
\r
1393 resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
\r
1394 thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
\r
1395 aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
\r
1396 resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
\r
1397 resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns
\r
1402 TGraph* graph = new TGraph ( kNptBins, fTransMomenta, resolution ) ;
\r
1405 graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
\r
1406 graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;
\r
1408 graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
\r
1409 graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;
\r
1411 graph->SetMinimum(1) ;
\r
1412 graph->SetMaximum(300.1) ;
\r
1413 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1414 graph->GetXaxis()->CenterTitle();
\r
1415 graph->GetXaxis()->SetNoExponent(1) ;
\r
1416 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1417 graph->GetYaxis()->CenterTitle();
\r
1419 graph->SetLineColor(color);
\r
1420 graph->SetMarkerColor(color);
\r
1421 graph->SetLineStyle(kDashed);
\r
1422 graph->SetLineWidth(linewidth);
\r
1428 TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
\r
1430 // particle = 0 ... choosen particle (setted particleMass)
\r
1431 // particle = 1 ... Pion
\r
1432 // particle = 2 ... Kaon
\r
1433 // particle = 3 ... D0
\r
1435 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1437 Double_t particleEfficiency[kNptBins]; // with chosen particle mass
\r
1438 Double_t kaonEfficiency[kNptBins], pionEfficiency[kNptBins], d0efficiency[kNptBins];
\r
1439 Double_t partEfficiency[2][400];
\r
1441 if (particle != 0) {
\r
1442 // resulting Pion and Kaon efficiency scaled with overall efficiency
\r
1443 Double_t doNotDecayFactor;
\r
1444 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
\r
1446 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1447 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
\r
1448 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
\r
1449 if ( massloop == 1 ) { // KAON
\r
1450 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
\r
1451 kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
\r
1453 doNotDecayFactor = 1.0 ;
\r
1454 pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
\r
1456 partEfficiency[0][j] = pionEfficiency[j];
\r
1457 partEfficiency[1][j] = kaonEfficiency[j];
\r
1461 // resulting estimate of the D0 efficiency
\r
1462 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1463 d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
\r
1466 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1467 particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
\r
1468 // NOTE: Decay factor (see kaon) should be included to be realiable
\r
1472 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1473 pionEfficiency[j] *= 100;
\r
1474 kaonEfficiency[j] *= 100;
\r
1475 d0efficiency[j] *= 100;
\r
1476 particleEfficiency[j] *= 100;
\r
1479 TGraph * graph = 0;
\r
1480 if (particle==0) {
\r
1481 graph = new TGraph ( kNptBins, fTransMomenta, particleEfficiency ) ; // choosen mass
\r
1482 graph->SetLineWidth(1);
\r
1483 } else if (particle==1) {
\r
1484 graph = new TGraph ( kNptBins, fTransMomenta, pionEfficiency ) ;
\r
1485 graph->SetLineWidth(1);
\r
1486 } else if (particle ==2) {
\r
1487 graph = new TGraph ( kNptBins, fTransMomenta, kaonEfficiency ) ;
\r
1488 graph->SetLineWidth(1);
\r
1489 } else if (particle ==3) {
\r
1490 graph = new TGraph ( kNptBins, fTransMomenta, d0efficiency ) ;
\r
1491 graph->SetLineStyle(kDashed);
\r
1495 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1496 graph->GetXaxis()->CenterTitle();
\r
1497 graph->GetXaxis()->SetNoExponent(1) ;
\r
1498 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1499 graph->GetYaxis()->SetTitle("Efficiency (%)") ;
\r
1500 graph->GetYaxis()->CenterTitle();
\r
1502 graph->SetMinimum(0.01) ;
\r
1503 graph->SetMaximum(100) ;
\r
1505 graph->SetLineColor(color);
\r
1506 graph->SetMarkerColor(color);
\r
1507 graph->SetLineWidth(linewidth);
\r
1512 TGraph * DetectorK::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {
\r
1514 // particle = 0 ... choosen particle (setted particleMass)
\r
1515 // particle = 1 ... Pion
\r
1516 // particle = 2 ... Kaon
\r
1519 Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1521 Double_t particleFake[kNptBins]; // with chosen particle mass
\r
1522 Double_t kaonFake[kNptBins], pionFake[kNptBins];
\r
1523 Double_t partFake[2][kNptBins];
\r
1525 if (particle != 0) {
\r
1526 // resulting Pion and Kaon efficiency scaled with overall efficiency
\r
1527 Double_t doNotDecayFactor;
\r
1528 for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
\r
1530 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1531 // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
\r
1532 Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
\r
1533 if ( massloop == 1 ) { // KAON
\r
1534 doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
\r
1535 kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;
\r
1537 pionFake[j] = fFake[0][j] ;
\r
1539 partFake[0][j] = pionFake[j];
\r
1540 partFake[1][j] = kaonFake[j];
\r
1545 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1546 particleFake[j] = fFake[2][j];
\r
1547 // NOTE: Decay factor (see kaon) should be included to be realiable
\r
1551 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1552 pionFake[j] *= 100;
\r
1553 kaonFake[j] *= 100;
\r
1554 particleFake[j] *= 100;
\r
1557 TGraph * graph = 0;
\r
1558 if (particle==0) {
\r
1559 graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
\r
1560 graph->SetLineWidth(1);
\r
1561 } else if (particle==1) {
\r
1562 graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
\r
1563 graph->SetLineWidth(1);
\r
1564 } else if (particle ==2) {
\r
1565 graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
\r
1566 graph->SetLineWidth(1);
\r
1569 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1570 graph->GetXaxis()->CenterTitle();
\r
1571 graph->GetXaxis()->SetNoExponent(1) ;
\r
1572 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1573 graph->GetYaxis()->SetTitle("Fake (%)") ;
\r
1574 graph->GetYaxis()->CenterTitle();
\r
1576 graph->SetMinimum(0.01) ;
\r
1577 graph->SetMaximum(100) ;
\r
1579 graph->SetLineColor(color);
\r
1580 graph->SetMarkerColor(color);
\r
1581 graph->SetLineWidth(linewidth);
\r
1585 TGraph * DetectorK::GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth) {
\r
1587 // particle = 0 ... choosen particle (setted particleMass)
\r
1588 // particle = 1 ... Pion
\r
1589 // particle = 2 ... Kaon
\r
1592 // Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
\r
1594 Double_t particleFake[kNptBins]; // with chosen particle mass
\r
1595 Double_t kaonFake[kNptBins], pionFake[kNptBins];
\r
1596 // Double_t partFake[2][kNptBins];
\r
1598 if (particle != 0) {
\r
1599 cout <<" not implemented"<<endl;
\r
1602 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1603 particleFake[j] = fFake[2][j];
\r
1604 // NOTE: Decay factor (see kaon) should be included to be realiable
\r
1609 for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
\r
1610 pionFake[j] = (1-pionFake[j])*100;
\r
1611 kaonFake[j] = (1-kaonFake[j])*100;
\r
1612 particleFake[j] = (1-particleFake[j])*100;
\r
1615 TGraph * graph = 0;
\r
1616 if (particle==0) {
\r
1617 graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
\r
1618 graph->SetLineWidth(1);
\r
1619 } else if (particle==1) {
\r
1620 graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
\r
1621 graph->SetLineWidth(1);
\r
1622 } else if (particle ==2) {
\r
1623 graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
\r
1624 graph->SetLineWidth(1);
\r
1627 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1628 graph->GetXaxis()->CenterTitle();
\r
1629 graph->GetXaxis()->SetNoExponent(1) ;
\r
1630 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1631 graph->GetYaxis()->SetTitle("Purity (%)") ;
\r
1632 graph->GetYaxis()->CenterTitle();
\r
1634 graph->SetMinimum(0.01) ;
\r
1635 graph->SetMaximum(100) ;
\r
1637 graph->SetLineColor(color);
\r
1638 graph->SetMarkerColor(color);
\r
1639 graph->SetLineWidth(linewidth);
\r
1645 TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
\r
1647 // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
\r
1648 // mode 0: impact parameter (convolution of pointing and vertex resolution)
\r
1649 // mode 1: pointing resolution
\r
1650 // mode 2: vtx resolution
\r
1653 TGraph *graph = new TGraph();
\r
1655 // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
\r
1656 TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
\r
1657 TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
\r
1659 TGraph *trackRes = GetGraphPointingResolution(axis,1);
\r
1660 Double_t *pt = trackRes->GetX();
\r
1661 Double_t *trRes = trackRes->GetY();
\r
1662 for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
\r
1663 Double_t vtxRes = 0;
\r
1665 vtxRes = vtxResRPhi.Eval(pt[ip]);
\r
1667 vtxRes = vtxResZ.Eval(pt[ip]);
\r
1670 graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
\r
1671 else if (mode ==1)
\r
1672 graph->SetPoint(ip,pt[ip],trRes[ip]);
\r
1674 graph->SetPoint(ip,pt[ip],vtxRes);
\r
1677 graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
\r
1678 graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;
\r
1680 graph->SetMinimum(1) ;
\r
1681 graph->SetMaximum(300.1) ;
\r
1682 graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
\r
1683 graph->GetXaxis()->CenterTitle();
\r
1684 graph->GetXaxis()->SetNoExponent(1) ;
\r
1685 graph->GetXaxis()->SetMoreLogLabels(1) ;
\r
1686 graph->GetYaxis()->CenterTitle();
\r
1688 graph->SetLineColor(color);
\r
1689 graph->SetMarkerColor(color);
\r
1690 graph->SetLineWidth(linewidth);
\r
1696 TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
\r
1698 // returns graph according to the number
\r
1702 return GetGraphPointingResolution(0,color, linewidth); // dr
\r
1704 return GetGraphPointingResolution(1,color, linewidth); // dz
\r
1706 return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
\r
1708 return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
\r
1710 return GetGraphMomentumResolution(color, linewidth); // pt resolution
\r
1712 return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle
\r
1714 return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion
\r
1716 return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon
\r
1718 return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0
\r
1720 return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle
\r
1722 return GetGraphRecoFakes(1, color, linewidth); // Fake pion
\r
1724 return GetGraphRecoFakes(2, color, linewidth); // Fake kaon
\r
1726 printf(" Error: chosen graph number not valid\n");
\r
1732 void DetectorK::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon) {
\r
1734 // All New configuration with X0 = 0.3 and resolution = 4 microns
\r
1736 AddLayer((char*)"bpipe",2.0,0.0022); // beam pipe
\r
1737 AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
\r
1739 // new ideal Pixel properties?
\r
1740 Double_t x0 = 0.0050;
\r
1741 Double_t resRPhi = 0.0006;
\r
1742 Double_t resZ = 0.0006;
\r
1750 AddLayer((char*)"ddd1", 2.2 , x0, resRPhi, resZ);
\r
1751 AddLayer((char*)"ddd2", 3.8 , x0, resRPhi, resZ);
\r
1752 AddLayer((char*)"ddd3", 6.8 , x0, resRPhi, resZ);
\r
1753 AddLayer((char*)"ddd4", 12.4 , x0, resRPhi, resZ);
\r
1754 AddLayer((char*)"ddd5", 23.5 , x0, resRPhi, resZ);
\r
1755 AddLayer((char*)"ddd6", 39.6 , x0, resRPhi, resZ);
\r
1756 AddLayer((char*)"ddd7", 43.0 , x0, resRPhi, resZ);
\r
1759 AddTPC(0.1,0.1); // TPC
\r
1763 void DetectorK::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) {
\r
1765 // Numbers taken from
\r
1766 // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
\r
1767 // number for misalingment: private communication with Andrea Dainese
\r
1769 AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe
\r
1770 AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
\r
1771 AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2
\r
1772 AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2
\r
1776 AddTPC(0.1,0.1); // TPC
\r
1778 // Adding the ITS - current configuration
\r
1780 if (AlignResiduals==0) {
\r
1782 AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130);
\r
1783 AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130);
\r
1784 AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025);
\r
1785 AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025);
\r
1786 AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830);
\r
1787 AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830);
\r
1789 } else if (AlignResiduals==1) {
\r
1791 // tracking errors ...
\r
1792 // (Additional systematic errors due to misalignments) ...
\r
1793 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
\r
1794 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
\r
1796 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
\r
1797 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1798 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
\r
1799 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1800 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
\r
1801 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1802 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
\r
1803 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1804 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
\r
1805 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1806 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
\r
1807 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1809 } else if (AlignResiduals==2) {
\r
1811 // tracking errors ... PLUS ... module misalignment
\r
1813 // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
\r
1814 // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
\r
1816 // the ITS modules are misalignment with small gaussian smearings with
\r
1817 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
\r
1819 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
\r
1820 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1821 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
\r
1822 TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
\r
1823 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
\r
1824 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1825 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
\r
1826 TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
\r
1827 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
\r
1828 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1829 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
\r
1830 TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
\r
1834 // the ITS modules are misalignment with small gaussian smearings with
\r
1835 // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
\r
1836 // unknown in Z ????
\r
1838 AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
\r
1839 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
\r
1840 AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
\r
1841 TMath::Sqrt(0.0130*0.0130+0.000*0.000));
\r
1842 AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
\r
1843 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
\r
1844 AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
\r
1845 TMath::Sqrt(0.0025*0.0025+0.000*0.000));
\r
1846 AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
\r
1847 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
\r
1848 AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
\r
1849 TMath::Sqrt(0.0830*0.0830+0.000*0.000));
\r
1857 void DetectorK::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
\r
1859 // Produces the standard performace plots
\r
1864 TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);
\r
1867 c1->cd(1); gPad->SetGridx(); gPad->SetGridy();
\r
1869 TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);
\r
1870 eff->SetTitle("Efficiencies");
\r
1872 if (!onlyPionEff) {
\r
1873 GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
\r
1874 GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
\r
1876 c1->cd(2); gPad->SetGridx(); gPad->SetGridy();
\r
1877 gPad->SetLogy(); gPad->SetLogx();
\r
1878 GetGraphMomentumResolution(color,linewidth)->Draw("AL");
\r
1880 c1->cd(3); gPad->SetGridx(); gPad->SetGridy();
\r
1882 GetGraphPointingResolution(0,color,linewidth)->Draw("AL");
\r
1884 c1->cd(4); gPad->SetGridx(); gPad->SetGridy();
\r
1886 GetGraphPointingResolution(1,color,linewidth)->Draw("AL");
\r
1890 TVirtualPad *c1 = gPad->GetMother();
\r
1893 GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");
\r
1894 if (!onlyPionEff) {
\r
1895 GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
\r
1896 GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
\r
1898 c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");
\r
1900 c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");
\r
1902 c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");
\r
1909 Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir) const
\r
1911 // Get local X of the track position estimated at the radius lab radius r.
\r
1912 // The track curvature is accounted exactly
\r
1914 // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
\r
1915 // 0 - take the intersection closest to the current track position
\r
1916 // >0 - go along the track (increasing fX)
\r
1917 // <0 - go backward (decreasing fX)
\r
1919 // special case of R=0
\r
1920 if (r<kAlmost0) {x=0; return kTRUE;}
\r
1922 const double* pars = tr->GetParameter();
\r
1923 const Double_t &fy=pars[0], &sn = pars[2];
\r
1925 double fx = tr->GetX();
\r
1926 double crv = tr->GetC(bz);
\r
1927 if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
\r
1928 if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
\r
1929 double det = (r-fx)*(r+fx);
\r
1930 if (det<0) return kFALSE; // does not reach raduis r
\r
1932 if (dir==0) return kTRUE;
\r
1933 det = TMath::Sqrt(det);
\r
1934 if (dir>0) { // along the track direction
\r
1935 if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
\r
1936 else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
\r
1938 else if(dir>0) { // agains track direction
\r
1939 if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
\r
1940 else if (fy>det) return kFALSE; // track is against Y axis
\r
1943 else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
\r
1944 double det = (r-fy)*(r+fy);
\r
1945 if (det<0) return kFALSE; // does not reach raduis r
\r
1946 det = TMath::Sqrt(det);
\r
1948 x = fx>0 ? det : -det; // choose the solution requiring the smalest step
\r
1951 else if (dir>0) { // along the track direction
\r
1952 if (fx > det) return kFALSE; // current point is in on the right from the circle
\r
1953 else if (fx <-det) x = -det; // on the left
\r
1954 else x = det; // within the circle
\r
1956 else { // against the track direction
\r
1957 if (fx <-det) return kFALSE;
\r
1958 else if (fx > det) x = det;
\r
1962 else { // general case of straight line
\r
1963 double cs = TMath::Sqrt((1-sn)*(1+sn));
\r
1964 double xsyc = fx*sn-fy*cs;
\r
1965 double det = (r-xsyc)*(r+xsyc);
\r
1966 if (det<0) return kFALSE; // does not reach raduis r
\r
1967 det = TMath::Sqrt(det);
\r
1968 double xcys = fx*cs+fy*sn;
\r
1970 if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
\r
1971 else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
\r
1972 if (t>=-det) t += -det; // take minimal step giving t>0
\r
1973 else return kFALSE; // both solutions have negative t
\r
1975 else { // go in increasing fx direction. (t+-det < 0)
\r
1976 if (t<det) t -= det; // take minimal step giving t<0
\r
1977 else return kFALSE; // both solutions have positive t
\r
1983 // get center of the track circle
\r
1984 double tR = 1./crv; // track radius (for the moment signed)
\r
1985 double cs = TMath::Sqrt((1-sn)*(1+sn));
\r
1986 double x0 = fx - sn*tR;
\r
1987 double y0 = fy + cs*tR;
\r
1988 double r0 = TMath::Sqrt(x0*x0+y0*y0);
\r
1989 // printf("Xc:%+e Yc:%+e\n",x0,y0);
\r
1991 if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
\r
1992 tR = TMath::Abs(tR);
\r
1993 double tR2r0 = tR/r0;
\r
1994 double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
\r
1995 double det = (1.-g)*(1.+g);
\r
1996 if (det<0) return kFALSE; // does not reach raduis r
\r
1997 det = TMath::Sqrt(det);
\r
1999 // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
\r
2000 // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
\r
2001 // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
\r
2003 double tmp = 1.+g*tR2r0;
\r
2005 double y = y0*tmp;
\r
2006 if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
\r
2007 double dfx = tR2r0*TMath::Abs(y0)*det;
\r
2008 double dfy = tR2r0*x0*TMath::Sign(det,y0);
\r
2009 if (dir==0) { // chose the one which corresponds to smallest step
\r
2010 double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
\r
2011 if (delta<0) x += dfx;
\r
2014 else if (dir>0) { // along track direction: x must be > fx
\r
2015 x -= dfx; // try the smallest step (dfx is positive)
\r
2016 if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;
\r
2018 else { // backward: x must be < fx
\r
2019 x += dfx; // try the smallest step (dfx is positive)
\r
2020 if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;
\r
2023 else { // special case: track touching the circle just in 1 point
\r
2024 if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE;
\r
2033 Double_t* DetectorK::PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay) {
\r
2036 printf("Error: Layer tracking efficiencies not set \n");
\r
2040 TMatrixD &tProbKomb = *probKomb;
\r
2041 TMatrixD &tProbLay = *probLay;
\r
2044 // Int_t base = tProbLay.GetNcols(); // 3? null, fake, correct
\r
2045 Int_t nLayer = tProbKomb.GetNcols(); // nlayer? - number of ITS layers
\r
2046 Int_t komb = tProbKomb.GetNrows(); // 3^nlayer? - number of kombinations
\r
2048 // Fill probabilities
\r
2050 Double_t probEff =0;
\r
2051 Double_t probFake =0;
\r
2052 for (Int_t num=0; num<komb; num++) {
\r
2053 Int_t flCorr=0, flFake=0, flNull=0;
\r
2054 for (Int_t l=0; l<nLayer; l++) {
\r
2055 if (tProbKomb(num,l)==0)
\r
2057 else if (tProbKomb(num,l)==1)
\r
2059 else if (tProbKomb(num,l)==2)
\r
2062 printf("Error: unexpected values in combinatorics table\n");
\r
2065 Int_t fkAtLeastCorr = fAtLeastCorr;
\r
2066 if (fAtLeastCorr == -1) fkAtLeastCorr = nLayer; // all hits are "correct"
\r
2068 if (flCorr>=fkAtLeastCorr && flFake==0) { // at least correct but zero fake
\r
2069 Double_t probEffLayer = 1;
\r
2070 for (Int_t l=0; l<nLayer; l++) {
\r
2071 probEffLayer *= tProbLay(tProbKomb(num,l),l);
\r
2072 // cout<<a(num,l)<<" ";
\r
2075 probEff+=probEffLayer;
\r
2078 if (flFake>=fAtLeastFake) {
\r
2079 Double_t probFakeLayer = 1;
\r
2080 for (Int_t l=0; l<nLayer; l++) {
\r
2081 probFakeLayer *= tProbLay(tProbKomb(num,l),l);
\r
2082 // cout<<a(num,l)<<" ";
\r
2085 probFake+=probFakeLayer;
\r
2089 Double_t *probs = new Double_t(2);
\r
2090 probs[0] = probEff; probs[1] = probFake;
\r