-#include "DetectorK.h"\r
-#include <TMath.h>\r
-#include <TMatrixD.h>\r
-#include <TGraph.h>\r
-#include <TAxis.h>\r
-#include <TFormula.h>\r
-#include <TCanvas.h>\r
-#include <TEllipse.h>\r
-#include <TText.h>\r
-#include <TGraphErrors.h>\r
-\r
-#include "AliExternalTrackParam.h"\r
-\r
-/***********************************************************\r
-\r
-Fast Simulation tool for Inner Tracker Systems\r
-\r
-original code of using the billoir technique was developed\r
-for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov\r
-http://rnc.lbl.gov/~jhthomas\r
-\r
-Changes by S. Rossegger -> see header file\r
-\r
-***********************************************************/\r
-\r
-\r
-#define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector\r
-\r
-#define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )\r
-#define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm)\r
-#define dNdEtaMinB 950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)\r
-#define dNdEtaCent 2300 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known)\r
-\r
-#define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)\r
-#define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors\r
-#define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)\r
-#define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0)\r
-#define EfficiencySearchFlag 2 // Define search method:\r
- // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0. \r
-\r
-#define PionMass 0.139 // Mass of the Pion\r
-#define KaonMass 0.498 // Mass of the Kaon\r
-#define D0Mass 1.865 // Mass of the D0\r
-\r
-\r
-\r
-class CylLayerK : public TNamed {\r
-public:\r
-\r
- CylLayerK(char *name) : TNamed(name,name) {}\r
- \r
- Float_t GetRadius() const {return radius;}\r
- Float_t GetRadL() const {return radL;}\r
- Float_t GetPhiRes() const {return phiRes;}\r
- Float_t GetZRes() const {return zRes;}\r
- Float_t GetLayerEff() const {return eff;}\r
-\r
- // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }\r
- Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes; \r
- Float_t eff;\r
- Bool_t isDead;\r
-\r
- ClassDef(CylLayerK,1);\r
-};\r
-\r
-\r
-class ForwardLayer : public TNamed {\r
-public:\r
- ForwardLayer(char *name) : TNamed(name,name) {}\r
- \r
- Float_t GetZ() const {return zPos;}\r
- Float_t GetXRes() const {return xRes;}\r
- Float_t GetYRes() const {return yRes;}\r
- Float_t GetThickness() const {return thickness;}\r
- Float_t Getdensity() const {return density;}\r
- Float_t GetLayerEff() const {return eff;}\r
-\r
- // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }\r
- Float_t zPos; Float_t xRes; Float_t yRes; \r
- Float_t radL;\r
- Float_t thickness;\r
- Float_t density;\r
- Float_t eff;\r
- Bool_t isDead;\r
-\r
- ClassDef(ForwardLayer,1);\r
-};\r
-\r
-\r
-ClassImp(DetectorK)\r
-DetectorK::DetectorK() \r
- : TNamed("test_detector","detector"),\r
- fNumberOfLayers(0),\r
- fNumberOfActiveLayers(0),\r
- fBField(0.5),\r
- fLhcUPCscale(1.0), \r
- fIntegrationTime(0.02), // in ms\r
- fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence\r
- fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle\r
- fParticleMass(0.140), // Standard: pion mass \r
- fMaxRadiusSlowDet(10.)\r
-{\r
- //\r
- // default constructor\r
- //\r
- // fLayers = new TObjArray();\r
- \r
-}\r
-\r
-DetectorK::DetectorK(char *name, char *title)\r
- : TNamed(name,title),\r
- fNumberOfLayers(0),\r
- fNumberOfActiveLayers(0),\r
- fBField(0.5),\r
- fLhcUPCscale(1.0),\r
- fIntegrationTime(0.02), // in ms\r
- fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence\r
- fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle\r
- fParticleMass(0.140), // Standard: pion mass\r
- fMaxRadiusSlowDet(10.)\r
-{\r
- //\r
- // default constructor, that set the name and title\r
- //\r
- // fLayers = new TObjArray();\r
-}\r
-DetectorK::~DetectorK() { // \r
- // virtual destructor\r
- //\r
- // delete fLayers;\r
-}\r
-\r
-void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) {\r
- //\r
- // Add additional layer to the list of layers (ordered by radius)\r
- // \r
-\r
- CylLayerK *newLayer = (CylLayerK*) fLayers.FindObject(name);\r
-\r
- if (!newLayer) {\r
- newLayer = new CylLayerK(name);\r
- newLayer->radius = radius;\r
- newLayer->radL = radL;\r
- newLayer->phiRes = phiRes;\r
- newLayer->zRes = zRes;\r
- newLayer->eff = eff;\r
-\r
- if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS) \r
- newLayer->isDead = kTRUE;\r
- else \r
- newLayer->isDead = kFALSE;\r
- \r
- if (fLayers.GetEntries()==0) \r
- fLayers.Add(newLayer);\r
- else {\r
- \r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- CylLayerK *l = (CylLayerK*)fLayers.At(i);\r
- if (radius<l->radius) {\r
- fLayers.AddBefore(l,newLayer);\r
- break;\r
- }\r
- if (radius>l->radius && (i+1)==fLayers.GetEntries() ) { \r
- // even bigger then last one\r
- fLayers.Add(newLayer);\r
- }\r
- }\r
- \r
- }\r
- fNumberOfLayers += 1;\r
- if (!(newLayer->isDead)) fNumberOfActiveLayers += 1;\r
-\r
-\r
- } else {\r
- printf("Layer with the name %s does already exist\n",name);\r
- }\r
- \r
-\r
-}\r
-\r
-void DetectorK::KillLayer(char *name) {\r
- //\r
- // Marks layer as dead. Contribution only by Material Budget\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot mark as dead\n",name);\r
- else {\r
- tmp->phiRes = 999999;\r
- tmp->zRes = 999999;\r
- if (!(tmp->isDead)) {\r
- tmp->isDead = kTRUE;\r
- fNumberOfActiveLayers -= 1; \r
- } \r
- }\r
-}\r
-\r
-void DetectorK::SetRadius(char *name, Float_t radius) {\r
- //\r
- // Set layer radius [cm]\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- \r
-\r
- if (!tmp) {\r
- printf("Layer %s not found - cannot set radius\n",name);\r
- } else {\r
- \r
- Float_t tmpRadL = tmp->radL;\r
- Float_t tmpPhiRes = tmp->phiRes;\r
- Float_t tmpZRes = tmp->zRes;\r
-\r
- RemoveLayer(name); // so that the ordering is correct\r
- AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes);\r
- }\r
-}\r
-\r
-Float_t DetectorK::GetRadius(char *name) {\r
- //\r
- // Return layer radius [cm]\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get radius\n",name);\r
- else \r
- return tmp->radius;\r
-\r
- return 0;\r
-}\r
-\r
-void DetectorK::SetRadiationLength(char *name, Float_t radL) {\r
- //\r
- // Set layer material [cm]\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set layer material\n",name);\r
- else {\r
- tmp->radL = radL;\r
- }\r
-}\r
-\r
-Float_t DetectorK::GetRadiationLength(char *name) {\r
- //\r
- // Return layer radius [cm]\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get layer material\n",name);\r
- else \r
- return tmp->radL;\r
- \r
- return 0;\r
- \r
-}\r
-\r
-void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) {\r
- //\r
- // Set layer resolution in [cm]\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set resolution\n",name);\r
- else {\r
-\r
- Bool_t wasDead = tmp->isDead;\r
- \r
- tmp->phiRes = phiRes;\r
- tmp->zRes = zRes;\r
- \r
- if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {\r
- tmp->isDead = kTRUE;\r
- if (!wasDead) fNumberOfActiveLayers -= 1;\r
- } else {\r
- tmp->isDead = kFALSE;\r
- if (wasDead) fNumberOfActiveLayers += 1;\r
- }\r
-\r
-\r
- }\r
-}\r
-\r
-Float_t DetectorK::GetResolution(char *name, Int_t axis) {\r
- //\r
- // Return layer resolution in [cm]\r
- // axis = 0: resolution in rphi\r
- // axis = 1: resolution in z\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get resolution\n",name);\r
- else {\r
- if (axis==0) return tmp->phiRes;\r
- if (axis==1) return tmp->zRes;\r
- printf("error: axis must be either 0 or 1 (rphi or z axis)\n");\r
- }\r
- return 0;\r
-}\r
-\r
-void DetectorK::SetLayerEfficiency(char *name, Float_t eff) {\r
- //\r
- // Set layer efficnecy (prop that his is missed within this layer) \r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set layer efficiency\n",name);\r
- else {\r
- tmp->eff = eff;\r
- }\r
-}\r
-\r
-Float_t DetectorK::GetLayerEfficiency(char *name) {\r
- //\r
- // Get layer efficnecy (prop that his is missed within this layer) \r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get layer efficneicy\n",name);\r
- else \r
- return tmp->eff;\r
- \r
- return 0;\r
- \r
-}\r
-\r
-void DetectorK::RemoveLayer(char *name) {\r
- //\r
- // Removes a layer from the list\r
- //\r
-\r
- CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot remove it\n",name);\r
- else {\r
- Bool_t wasDead = tmp->isDead;\r
- fLayers.Remove(tmp);\r
- fNumberOfLayers -= 1;\r
- if (!wasDead) fNumberOfActiveLayers -= 1;\r
- }\r
-}\r
-\r
-\r
-void DetectorK::PrintLayout() {\r
- //\r
- // Prints the detector layout\r
- //\r
-\r
- printf("Detector %s: \"%s\"\n",GetName(),GetTitle());\r
- \r
- if (fLayers.GetEntries()>0) \r
- printf(" Name \t\t r [cm] \t X0 \t phi & z res [um]\n");\r
-\r
- CylLayerK *tmp = 0;\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- tmp = (CylLayerK*)fLayers.At(i);\r
- \r
- // don't print all the tpc layers\r
- TString name(tmp->GetName());\r
- if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
-\r
- printf("%d. %s \t %03.2f \t%1.4f\t ",i,\r
- tmp->GetName(), tmp->radius, tmp->radL);\r
- if (tmp->phiRes==RIDICULOUS) \r
- printf(" - ");\r
- else\r
- printf("%3.0f ",tmp->phiRes*10000);\r
- if (tmp->zRes==RIDICULOUS) \r
- printf(" -\n");\r
- else\r
- printf("%3.0f\n",tmp->zRes*10000);\r
- }\r
-}\r
-\r
-void DetectorK::PlotLayout(Int_t plotDead) {\r
- //\r
- // Plots the detector layout in Front view\r
- //\r
-\r
- Double_t x0=0, y0=0;\r
-\r
- TGraphErrors *gr = new TGraphErrors();\r
- gr->SetPoint(0,0,0);\r
- CylLayerK *lastLayer = (CylLayerK*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->radius;\r
- gr->SetPointError(0,maxRad,maxRad);\r
- gr->Draw("APE");\r
- \r
-\r
- CylLayerK *tmp = 0;\r
- for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {\r
- tmp = (CylLayerK*)fLayers.At(i);\r
- \r
-\r
- Double_t txtpos = tmp->radius;\r
- if ((tmp->isDead)) txtpos*=-1; //\r
- TText *txt = new TText(x0,txtpos,tmp->GetName());\r
- txt->SetTextSizePixels(5); txt->SetTextAlign(21);\r
- if (!tmp->isDead || plotDead) txt->Draw();\r
-\r
- TEllipse *layEl = new TEllipse(x0,y0,tmp->radius);\r
- // layEl->SetFillColor(5);\r
- layEl->SetFillStyle(5001);\r
- layEl->SetLineStyle(tmp->isDead+1); // dashed if not active\r
- layEl->SetLineColor(4);\r
- TString name(tmp->GetName());\r
- if (!tmp->isDead) layEl->SetLineWidth(2);\r
- if (name.Contains("tpc") ) layEl->SetLineColor(29);\r
-\r
- if (!tmp->isDead || plotDead) layEl->Draw();\r
- \r
- }\r
-\r
-}\r
-\r
-\r
-\r
-void DetectorK::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {\r
- //\r
- // Emulates the TPC\r
- // \r
- // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow \r
-\r
-\r
- AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage\r
- \r
- // % Radiation Lengths ... Average per TPC row (i.e. total/159 )\r
- Float_t radLPerRow = 0.000036;\r
- \r
- Float_t tpcInnerRadialPitch = 0.75 ; // cm\r
- Float_t tpcMiddleRadialPitch = 1.0 ; // cm\r
- Float_t tpcOuterRadialPitch = 1.5 ; // cm\r
- // Float_t tpcInnerPadWidth = 0.4 ; // cm\r
- // Float_t tpcMiddlePadWidth = 0.6 ; // cm\r
- // Float_t tpcOuterPadWidth = 0.6 ; // cm\r
- Float_t innerRows = 63 ;\r
- Float_t middleRows = 64 ;\r
- Float_t outerRows = 32 ;\r
- Float_t tpcRows = (innerRows + middleRows + outerRows) ;\r
- Float_t rowOneRadius = 85.2 ; // cm\r
- Float_t row64Radius = 135.1 ; // cm\r
- Float_t row128Radius = 199.2 ; // cm \r
- \r
- for ( Int_t k = 0 ; k < tpcRows ; k++ ) {\r
- \r
- Float_t rowRadius =0;\r
- if (k<innerRows) \r
- rowRadius = rowOneRadius + k*tpcInnerRadialPitch ;\r
- else if ( k>=innerRows && k<(innerRows+middleRows) )\r
- rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;\r
- else if (k>=(innerRows+middleRows) && k<tpcRows )\r
- rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;\r
-\r
- if ( k%skip == 0 )\r
- AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean); \r
- else \r
- AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row\r
- \r
- \r
- }\r
- \r
-}\r
-\r
-void DetectorK::RemoveTPC() {\r
-\r
- // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)\r
- CylLayerK *tmp = 0;\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- tmp = (CylLayerK*)fLayers.At(i); \r
- TString name(tmp->GetName());\r
- if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }\r
- }\r
- RemoveLayer((char*)"IFC");\r
- \r
-}\r
-\r
-\r
-Double_t DetectorK::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const\r
-{\r
- //\r
- // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)\r
- //\r
-\r
- Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ;\r
- Double_t theta = 0.0 ; // Momentum and mass in GeV\r
- // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );\r
- if ( radLength > 0 ) theta = 0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ;\r
- return (theta) ;\r
-}\r
-\r
-\r
-Double_t DetectorK::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
-{\r
- // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm \r
- // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm\r
- // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius\r
- Double_t sx, sy, goodHit ;\r
- sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;\r
- sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ;\r
- goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ;\r
- return ( goodHit ) ;\r
-}\r
-\r
-\r
-Double_t DetectorK::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
-{\r
- // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm\r
- // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function\r
- Double_t sx, goodHit ;\r
- sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;\r
- goodHit = 1./(1+sx) ;\r
- return ( goodHit ) ; \r
-}\r
-\r
-Double_t DetectorK::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ ) \r
-{\r
- // Based on work by Ruben Shahoyen \r
- // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function\r
- // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account \r
- // Following is correct for 2 DOF\r
-\r
- Double_t gamma = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level\r
- Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2; \r
- Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*gamma));\r
-\r
- return ( goodHit ) ; \r
-}\r
-\r
-\r
-Double_t DetectorK::HitDensity ( Double_t radius ) \r
-{\r
- // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.\r
- // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results\r
- // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.\r
- // Note that this function assumes we are working in CM and CM**2 [not meters].\r
- // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.\r
-\r
- // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors \r
- // and only fast detectors reside outside this radius.\r
- Double_t arealDensity = 0 ;\r
-\r
- if ( radius > fMaxRadiusSlowDet ) \r
- {\r
- arealDensity = OneEventHitDensity(dNdEtaCent,radius) ; // Fast detectors see central collision density (only)\r
- arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background \r
- }\r
-\r
- if (radius < fMaxRadiusSlowDet )\r
- { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.\r
- arealDensity = OneEventHitDensity(dNdEtaCent,radius) \r
- + IntegratedHitDensity(dNdEtaMinB,radius) \r
- + UpcHitDensity(radius) ;\r
- arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ; \r
- // Increase density due to background \r
- } \r
-\r
- return ( arealDensity ) ; \r
-}\r
-\r
-\r
-double DetectorK::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const\r
-{\r
- // This is for one event at the vertex. No smearing.\r
- double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?\r
- // note: surface of sphere is '4*pi*r^2'\r
- // surface of cylinder is '2*pi*r* h' \r
- return den ;\r
-} \r
-\r
-\r
-double DetectorK::IntegratedHitDensity(Double_t multiplicity, Double_t radius)\r
-{ \r
- // The integral of minBias events smeared over a gaussian vertex distribution.\r
- // Based on work by Yan Lu 12/20/2006, all radii in centimeters.\r
-\r
- Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;\r
- Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;\r
-\r
- // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.\r
- if ( den < OneEventHitDensity(multiplicity,radius) ) den = OneEventHitDensity(multiplicity,radius) ; \r
-\r
- return den ;\r
-} \r
-\r
-\r
-double DetectorK::UpcHitDensity(Double_t radius)\r
-{ \r
- // QED electrons ...\r
-\r
- Double_t mUPCelectrons ; ; \r
- // mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC\r
- mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB; // Fit to 'Rossegger,Sadovsky'-Alice simulation\r
- if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R\r
- mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate\r
- mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC\r
-\r
- return mUPCelectrons ;\r
-} \r
-\r
-\r
-double DetectorK::Dist(double z, double r)\r
-{\r
- // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution\r
- // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm\r
- // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.\r
- Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration\r
- Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work\r
- double dist = 0.0 ;\r
- double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm\r
- double z0 = 0.0 ; //cm\r
- for(int i=0; i<nsteps; i++){\r
- if ( i == nsteps-1 ) index = 1 ;\r
- z0 = -4*SigmaD + i*dz0 ;\r
- dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) * \r
- (1/sqrt((z-z0)*(z-z0) + r*r)) ;\r
- if ( index != 4 ) index = 4; else index = 2 ;\r
- }\r
- return dist; \r
-}\r
-\r
-#define PZero 0.861 // Momentum of back to back decay particles in the CM frame\r
-#define EPiZero 0.872 // Energy of the pion from a D0 decay at rest\r
-#define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest\r
-\r
-Double_t DetectorK::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const {\r
- // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi.\r
-\r
- Double_t const1 = pt / D0Mass ;\r
- Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;\r
- Double_t sum, ptPi, ptK ;\r
- Double_t effp, effk ;\r
-\r
- sum = 0.0 ;\r
- for ( Int_t k = 0 ; k < 360 ; k++ ) {\r
- \r
- Double_t theta = k * TMath::Pi() / 180. ;\r
- \r
- ptPi = TMath::Sqrt( \r
- PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +\r
- const1*const1*EPiZero*EPiZero -\r
- 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +\r
- PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)\r
- ) ;\r
- \r
- ptK = TMath::Sqrt( \r
- PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +\r
- const1*const1*EKZero*EKZero +\r
- 2*PZero*TMath::Cos(theta)*const2*const1*EKZero +\r
- PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)\r
- ) ;\r
-\r
- // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays\r
- Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; \r
- Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; \r
- \r
- if ( pionindex >= 400 ) pionindex = 399 ;\r
- if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ;\r
- if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd\r
- if ( effp < 0 ) effp = 0 ;\r
-\r
- if ( kaonindex >= 400 ) kaonindex = 399 ;\r
- if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ;\r
- if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd\r
- if ( effk < 0 ) effk = 0 ;\r
-\r
- // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.\r
- \r
- sum += effp * effk ;\r
- \r
- } \r
- \r
- Double_t mean =sum/360; \r
- return mean ;\r
- \r
-}\r
-\r
-\r
-void DetectorK::SolveDOFminusOneAverage() {\r
- // \r
- // Short study to address "# layers-1 efficiencies"\r
- // saves the means in the according arrays\r
- // Note: Obviously, does not work for the Telescope equation \r
- //\r
- \r
- Double_t fMomentumResM[400], fResolutionRPhiM[400], fResolutionZM[400]; \r
- Double_t efficiencyM[3][400];\r
- for (Int_t i=0; i<400; i++) {\r
- fMomentumResM[i] = 0; // Momentum resolution\r
- fResolutionRPhiM[i] = 0; // Resolution in R\r
- fResolutionZM[i] = 0; // Resolution in Z\r
- for (Int_t part=0; part<3; part++) \r
- efficiencyM[part][i] = 0; // efficiencies\r
- }\r
-\r
- // loop over active layers in ITS (remove 1 by 1)\r
- Int_t nITSLayers = 0;\r
- CylLayerK *layer =0;\r
- for (Int_t j=0; j<(fLayers.GetEntries()-1); j++) { \r
- layer = (CylLayerK*)fLayers.At(j);\r
- TString name(layer->GetName());\r
- if (name.Contains("tpc")) continue;\r
- if (!(layer->isDead)) {\r
-\r
- nITSLayers++; \r
- printf("Kill Layer %s\n",name.Data());\r
- Double_t rRes = GetResolution((char*)name.Data(),0);\r
- Double_t zRes = GetResolution((char*)name.Data(),1);\r
- KillLayer((char*)name.Data());\r
- // PrintLayout();\r
- SolveViaBilloir(1,0); \r
-\r
- // produce sum for the mean calculation\r
- for (Int_t i=0; i<400; i++) {\r
- fMomentumResM[i] += fMomentumRes[i]; // Momentum resolution\r
- fResolutionRPhiM[i] += fResolutionRPhi[i]; // Resolution in R\r
- fResolutionZM[i] += fResolutionZ[i]; // Resolution in Z\r
- for (Int_t part=0; part<3; part++) \r
- efficiencyM[part][i] += fEfficiency[part][i]; // efficiencies\r
- }\r
-\r
- // "Restore" layer ...\r
- SetResolution((char*)name.Data(),rRes,zRes); \r
- \r
- }\r
- }\r
- \r
- // save means in "std. Arrays"\r
- for (Int_t i=0; i<400; i++) {\r
- fMomentumRes[i] = fMomentumResM[i]/nITSLayers; // Momentum resolution\r
- fResolutionRPhi[i] = fResolutionRPhiM[i]/nITSLayers; // Resolution in R\r
- fResolutionZ[i] = fResolutionZM[i]/nITSLayers; // Resolution in Z\r
- for (Int_t part=0; part<3; part++) \r
- fEfficiency[part][i] = efficiencyM[part][i]/nITSLayers; // efficiencies\r
- }\r
-\r
-\r
-}\r
-\r
-void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt) {\r
- //\r
- // Solves the current geometry with the Billoir technique \r
- // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. )\r
- // ABOVE IS OBSOLETE -> NOW, its uses the Aliroot Kalman technique\r
- //\r
-\r
- static AliExternalTrackParam probTr; // track to propagate\r
- probTr.SetUseLogTermMS(kTRUE);\r
-\r
-\r
- Int_t nPt = 400;\r
- // Clean up ......\r
- for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {\r
- for (Int_t j=0; j<nPt; j++) {\r
- fDetPointRes[i][j] = RIDICULOUS;\r
- fDetPointZRes[i][j] = RIDICULOUS;\r
- fTransMomenta[i] =0;\r
- fMomentumRes[i] =0;\r
- fResolutionRPhi[i] =0;\r
- }\r
- }\r
- \r
- if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt'\r
- nPt = 3;\r
- }\r
-\r
-\r
- // Calculate track parameters using Billoirs method of matrices\r
-\r
- Double_t pt,tgl, pz, lambda, deltaPoverP ;\r
- Double_t charge ;\r
- Double_t mass[3] ;\r
- Int_t printOnce = 1 ;\r
-\r
- mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k \r
-\r
- mass[2] = fParticleMass; // third loop\r
-\r
- Int_t mStart =0; \r
- if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode\r
-\r
- for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) { \r
- \r
- // PseudoRapidity OK, used as an angle\r
- lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ; \r
- \r
-\r
- for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop\r
- \r
- CylLayerK *last = (CylLayerK*) fLayers.At((fLayers.GetEntries()-1));\r
-\r
- // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV\r
- Double_t bigRad = last->radius/2 ; // min. pt which the algorithm below could handle\r
- //if (bigRad<61) bigRad=61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) )\r
- fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 + TMath::Power(10,2.3*i/nPt) / 10.0 ; \r
- if (!allPt) { // just 3 points around meanPt\r
- fTransMomenta[i] = meanPt-0.1+i*0.1;\r
- }\r
- \r
-\r
- // New from here ................\r
-\r
- // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis\r
- // These are the EndPoint values for y, z, a, b, and d\r
- double bGauss = fBField*10; // field in kgauss\r
- pt = fTransMomenta[i]; // GeV/c\r
- tgl = TMath::Tan(lambda); // dip\r
- charge = -1; // Assume an electron \r
- pz = pt * TMath::Tan(lambda) ; // GeV/\r
- enum {kY,kZ,kSnp,kTgl,kPtI}; // track parameter aliases\r
- enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases\r
- //\r
- probTr.Reset();\r
- double *trPars = (double*)probTr.GetParameter();\r
- double *trCov = (double*)probTr.GetCovariance();\r
- trPars[kY] = 0; // start from Y = 0\r
- trPars[kZ] = 0; // Z = 0 \r
- trPars[kSnp] = 0; // track along X axis at the vertex\r
- trPars[kTgl] = TMath::Tan(lambda); // dip\r
- trPars[kPtI] = charge/pt; // q/pt\r
- //\r
- // put tiny errors to propagate to the outer radius\r
- trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;\r
- double xR = 0;\r
- if (!GetXatLabR(&probTr, last->radius ,xR,bGauss,1)) {\r
- printf("Track with pt=%f cannot reach radius %f\n",pt,last->radius);\r
- continue;\r
- }\r
- probTr.PropagateTo(xR, bGauss); // bring track to outer layer\r
- // reset cov.matrix\r
- const double kLargeErr2Coord = 50*50;\r
- const double kLargeErr2Dir = 0.6*0.6;\r
- const double kLargeErr2PtI = 0.5*0.5;\r
- for (int ic=15;ic--;) trCov[ic] = 0.;\r
- trCov[kY2] = trCov[kZ2] = kLargeErr2Coord; \r
- trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;\r
- trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];\r
- //\r
- // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);\r
-\r
- // Set Detector-Efficiency Storage area to unity\r
- fEfficiency[massloop][i] = 1.0 ;\r
- //\r
- // Back-propagate the covariance matrix along the track. \r
- \r
- CylLayerK *layer = 0;\r
- \r
- // find last "active layer" - start tracking at the last active layer \r
- Int_t lastActiveLayer = 0;\r
- for (Int_t j=fLayers.GetEntries(); j--;) { \r
- layer = (CylLayerK*)fLayers.At(j);\r
- if (!(layer->isDead)) { // is alive\r
- lastActiveLayer = j;\r
- break;\r
- }\r
- }\r
- \r
- for (Int_t j=lastActiveLayer; j--;) { // Layer loop\r
-\r
- layer = (CylLayerK*)fLayers.At(j);\r
- TString name(layer->GetName());\r
- Bool_t isVertex = name.Contains("vertex");\r
- //\r
- if (!GetXatLabR(&probTr, layer->radius ,xR,bGauss)) { \r
- printf("Track with pt=%f cannot reach radius %f. This should not happen here\n",pt,layer->radius);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- probTr.PropagateTo(xR, bGauss); // propagate to this layer\r
- //\r
- // rotate to frame with X axis normal to the surface\r
- if (!isVertex) {\r
- double pos[3];\r
- probTr.GetXYZ(pos); // lab position\r
- double phi = TMath::ATan2(pos[1],pos[0]);\r
- if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);\r
- if (!probTr.Rotate(phi)) {\r
- printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f\n",\r
- phi,layer->radius,pos[0],pos[1],pos[2]);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- }\r
- // save resolutions at this layer\r
- fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters\r
- fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters\r
- //printf(">> L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);\r
- // End save\r
- //\r
- if (isVertex) continue;\r
- //\r
- // create fake measurement with the errors assigned to the layer\r
- // account for the measurement there \r
- double meas[2] = {probTr.GetY(),probTr.GetZ()};\r
- double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};\r
- //\r
- if (!probTr.Update(meas,measErr2)) {\r
- printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",\r
- meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);\r
- probTr.Print();\r
- exit(1);\r
- }\r
-\r
- // correct for materials of this layer\r
- // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param\r
- if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {\r
- printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- //\r
- }\r
- \r
- // Pattern recognition is done .... save values like vertex resolution etc.\r
-\r
- // Convert the Convariance matrix parameters into physical quantities\r
- // The results are propogated to the previous point but *do not* include the measurement at that point.\r
- // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P(); // Absolute magnitude so ignore charge\r
- deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P(); \r
- fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent\r
- fResolutionRPhi[i] = TMath::Sqrt( probTr.GetSigmaY2() ) * 1.e4; // result in microns\r
- fResolutionZ[i] = TMath::Sqrt( probTr.GetSigmaZ2() ) * 1.e4; // result in microns\r
- // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius\r
- \r
- \r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
- printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;\r
- printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);\r
- printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ;\r
- // printOnce =0;\r
- }\r
- \r
- // print out and efficiency calculation\r
- // for (Int_t j=fLayers.GetEntries(); j--;) { // Layer loop\r
- for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop\r
- \r
- layer = (CylLayerK*)fLayers.At(j);\r
- \r
- // Convert to Meters, Tesla, and GeV\r
- Float_t radius = layer->radius /100;\r
- Float_t phiRes = layer->phiRes /100;\r
- Float_t zRes = layer->zRes /100;\r
- Float_t radLength = layer->radL;\r
- Float_t leff = layer->eff; // basic layer efficiency\r
- Bool_t isDead = layer->isDead;\r
- \r
-\r
- if ( (!isDead && radLength >0) ) { \r
- Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + \r
- phiRes * phiRes ) * 100. ; // work in cm\r
- Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +\r
- zRes * zRes ) * 100. ; // work in cm\r
- \r
- Double_t layerEfficiency = 0;\r
- if ( EfficiencySearchFlag == 0 )\r
- layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;\r
- else if ( EfficiencySearchFlag == 1 )\r
- layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;\r
- else if ( EfficiencySearchFlag == 2 )\r
- layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;\r
-\r
- TString name(layer->GetName());\r
- if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
-\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
-\r
-\r
- printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",\r
- layer->GetName(), radius*100, radLength, \r
- fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,\r
- phiRes*1.e6, zRes*1.e6,\r
- HitDensity(radius*100)) ;\r
- if (!name.Contains("tpc")) \r
- printf("%10.3f\n", layerEfficiency);\r
- else\r
- printf(" - \n");\r
- }\r
-\r
- if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;\r
-\r
- }\r
- /*\r
- // vertex print\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {\r
- printf("%s:\t ----- ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6);\r
- }\r
- */\r
- }\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
- if (fNumberOfActiveLayers >=15) printOnce = 0 ;\r
- printf("\n") ;\r
- }\r
- \r
-\r
-\r
-\r
- if (fNumberOfActiveLayers <15 ) {\r
-\r
-\r
-\r
- // BACKWORD TRACKING +++++++++++++++++\r
- // number of layers is quite low ... efficiency calculation was probably nonsense \r
- // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates"\r
-\r
- // For below, see paper, NIM A262 (1987) p.444, eqs.12.\r
- // Equivalently, one can simply combine the forward and backward estimates. Assuming\r
- // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can\r
- // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1).\r
- // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.\r
-\r
- Bool_t doLikeAliRoot = kFALSE; // don't do the "combined info" but do like in Aliroot\r
-\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) { \r
- printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n",\r
- fNumberOfActiveLayers); \r
- }\r
- \r
- // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam)\r
- // mIstar.UnitMatrix(); // start with unity\r
- if (doLikeAliRoot) {\r
- probTr.ResetCovariance(10.);\r
- } else {\r
- // cannot do complete reset, set to very large errors\r
- trCov[kY2] = trCov[kZ2] = kLargeErr2Coord; \r
- trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;\r
- trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];\r
- // cout<<pt<<": "<<kLargeErr2Coord<<" "<<kLargeErr2Dir<<" "<<kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]<<endl;\r
- }\r
- // Clean up and storing of "forward estimates"\r
- Double_t detPointResForw[kMaxNumberOfDetectors][400], detPointZResForw[kMaxNumberOfDetectors][400] ; \r
- for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {\r
- for (Int_t l=0; l<nPt; l++) {\r
- detPointResForw[k][l] = fDetPointRes[k][l];\r
- if (!doLikeAliRoot) fDetPointRes[k][l] = RIDICULOUS;\r
- detPointZResForw[k][l] = fDetPointZRes[k][l];\r
- if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS;\r
- }\r
- }\r
- \r
- // find first "active layer" - start tracking at the first active layer \r
- Int_t firstActiveLayer = 0;\r
- for (Int_t j=0; j<fLayers.GetEntries(); j++) { \r
- layer = (CylLayerK*)fLayers.At(j);\r
- if (!(layer->isDead)) { // is alive\r
- firstActiveLayer = j;\r
- break;\r
- }\r
- }\r
- probTr.Rotate(0);\r
- for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()); j++) { // Layer loop\r
- \r
- layer = (CylLayerK*)fLayers.At(j);\r
- // CylLayerK *nextlayer = (CylLayerK*)fLayers.At(j+1);\r
-\r
- TString name(layer->GetName());\r
- Bool_t isVertex = name.Contains("vertex");\r
- //\r
- if (!GetXatLabR(&probTr, layer->radius ,xR,bGauss)) { \r
- printf("Track with pt=%f cannot reach radius %f. This should not happen here\n",pt,layer->radius);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- probTr.PropagateTo(xR, bGauss); // propagate to this layer\r
- if (!isVertex) {\r
- // rotate to frame with X axis normal to the surface\r
- double pos[3];\r
- probTr.GetXYZ(pos); // lab position\r
- double phi = TMath::ATan2(pos[1],pos[0]);\r
- if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);\r
- if (!probTr.Rotate(phi)) {\r
- printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f\n",\r
- phi,layer->radius,pos[0],pos[1],pos[2]);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- }\r
- // \r
- fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters\r
- fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters\r
- //\r
- //printf("<< L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);\r
- // create fake measurement with the errors assigned to the layer\r
- // account for the measurement there\r
- if (isVertex) continue;\r
- double meas[2] = {probTr.GetY(),probTr.GetZ()};\r
- double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};\r
- //\r
- if (!probTr.Update(meas,measErr2)) {\r
- printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",\r
- meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- // correct for materials of this layer\r
- // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param\r
- if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {\r
- printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);\r
- probTr.Print();\r
- exit(1);\r
- }\r
- }\r
-\r
- // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!!\r
- // ++++++++++++++\r
- // also update the values for the track position ??????\r
- /*\r
- // Pattern recognition is done .... save values like vertex resolution etc.\r
- \r
- // Invert the Matrix to recover the convariance matrix\r
- mIstar.Invert() ;\r
- // Convert the Convariance matrix parameters into physical quantities\r
- // The results are propogated to the previous point but *do not* include the measurement at that point.\r
- deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge\r
- fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent\r
- fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns\r
- fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns\r
- // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius\r
- */\r
- \r
- // Weighted combination of the forward and backward estimates\r
- if (!doLikeAliRoot) {\r
- for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { \r
- fDetPointRes[j][i] = 1/(1/detPointResForw[j][i] + 1/fDetPointRes[j][i]); \r
- fDetPointZRes[j][i] = 1/(1/detPointZResForw[j][i] + 1/fDetPointZRes[j][i]); \r
- }\r
- }\r
- // Set Detector-Efficiency Storage area to unity\r
- fEfficiency[massloop][i] = 1.0 ;\r
- \r
- // print out and efficiency calculation\r
- for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop\r
- \r
- layer = (CylLayerK*)fLayers.At(j);\r
- \r
- // Convert to Meters, Tesla, and GeV\r
- Float_t radius = layer->radius /100;\r
- Float_t phiRes = layer->phiRes /100;\r
- Float_t zRes = layer->zRes /100;\r
- Float_t radLength = layer->radL;\r
- Float_t leff = layer->eff;\r
- Bool_t isDead = layer->isDead;\r
- \r
-\r
- if ( (!isDead && radLength >0) ) { \r
- Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] + \r
- phiRes * phiRes ) * 100. ; // work in cm\r
- Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +\r
- zRes * zRes ) * 100. ; // work in cm\r
- \r
- Double_t layerEfficiency = 0;\r
- if ( EfficiencySearchFlag == 0 )\r
- layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;\r
- else if ( EfficiencySearchFlag == 1 )\r
- layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;\r
- else if ( EfficiencySearchFlag == 2 )\r
- layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;\r
-\r
- TString name(layer->GetName());\r
- if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;\r
-\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
-\r
-\r
- printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",\r
- layer->GetName(), radius*100, radLength, \r
- fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,\r
- phiRes*1.e6, zRes*1.e6,\r
- HitDensity(radius*100)) ;\r
- if (!name.Contains("tpc")) \r
- printf("%10.3f\n", layerEfficiency);\r
- else\r
- printf(" - \n");\r
- }\r
-\r
- if (!name.Contains("tpc")) fEfficiency[massloop][i] *= layerEfficiency;\r
- \r
- }\r
- }\r
- if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {\r
- printOnce = 0 ;\r
- printf("\n") ;\r
- }\r
- } \r
- } // mass loop\r
- } // pt loop\r
-\r
- probTr.SetUseLogTermMS(kFALSE); // Reset of MS term usage to avoid problems since its static\r
-\r
-\r
-}\r
-\r
-\r
-TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {\r
- //\r
- // returns the momentum resolution \r
- //\r
- \r
- TGraph *graph = new TGraph(400, fTransMomenta, fMomentumRes);\r
- graph->SetTitle("Momentum Resolution .vs. Pt" ) ;\r
- // graph->GetXaxis()->SetRangeUser(0.,5.0) ;\r
- graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
- graph->GetXaxis()->CenterTitle();\r
- graph->GetXaxis()->SetNoExponent(1) ;\r
- graph->GetXaxis()->SetMoreLogLabels(1) ;\r
- graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;\r
- graph->GetYaxis()->CenterTitle();\r
-\r
- graph->SetMaximum(20) ;\r
- graph->SetMinimum(0.1) ;\r
- graph->SetLineColor(color);\r
- graph->SetMarkerColor(color);\r
- graph->SetLineWidth(linewidth);\r
-\r
- return graph;\r
-\r
-}\r
-\r
-TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {\r
- \r
- // Returns the pointing resolution\r
- // axis = 0 ... rphi pointing resolution\r
- // axis = 1 ... z pointing resolution\r
- //\r
-\r
- TGraph * graph = 0;\r
-\r
- if (axis==0) {\r
- graph = new TGraph ( 400, fTransMomenta, fResolutionRPhi ) ;\r
- graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;\r
- graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;\r
- } else {\r
- graph = new TGraph ( 400, fTransMomenta, fResolutionZ ) ;\r
- graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;\r
- graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;\r
- }\r
- \r
- graph->SetMinimum(1) ;\r
- graph->SetMaximum(300.1) ;\r
- graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
- graph->GetXaxis()->CenterTitle();\r
- graph->GetXaxis()->SetNoExponent(1) ;\r
- graph->GetXaxis()->SetMoreLogLabels(1) ;\r
- graph->GetYaxis()->CenterTitle();\r
- \r
- graph->SetLineWidth(linewidth);\r
- graph->SetLineColor(color);\r
- graph->SetMarkerColor(color);\r
- \r
- return graph;\r
-\r
-}\r
-\r
-\r
-TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {\r
- //\r
- // returns the Pointing resolution (accoring to Telescope equation)\r
- // axis =0 ... in rphi\r
- // axis =1 ... in z\r
- //\r
- \r
- Double_t resolution[400];\r
-\r
- Double_t layerResolution[2];\r
- Double_t layerRadius[2];\r
- Double_t layerThickness[2];\r
-\r
- Int_t count =0; // search two first active layers\r
- printf("Telescope equation for layers: ");\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- CylLayerK *l = (CylLayerK*)fLayers.At(i);\r
- if (!l->isDead && l->radius>0) {\r
- layerRadius[count] = l->radius;\r
- layerThickness[count] = l->radL;\r
- if (axis==0) {\r
- layerResolution[count] = l->phiRes;\r
- } else {\r
- layerResolution[count] = l->zRes;\r
- }\r
- printf("%s, ",l->GetName());\r
- count++;\r
- }\r
- if (count>=2) break; \r
- }\r
- printf("\n");\r
-\r
- Double_t pt, momentum, thickness,aMCS ;\r
- Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
-\r
- for ( Int_t i = 0 ; i < 400 ; i++ ) { \r
- // Reference data as if first two layers were acting all alone \r
- pt = fTransMomenta[i] ;\r
- momentum = pt / TMath::Cos(lambda) ; // Total momentum\r
- resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1] \r
- + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;\r
- resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;\r
- thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;\r
- aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;\r
- resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;\r
- resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns\r
- }\r
-\r
-\r
-\r
- TGraph* graph = new TGraph ( 400, fTransMomenta, resolution ) ;\r
- \r
- if (axis==0) {\r
- graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;\r
- graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;\r
- } else {\r
- graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;\r
- graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;\r
- }\r
- graph->SetMinimum(1) ;\r
- graph->SetMaximum(300.1) ;\r
- graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
- graph->GetXaxis()->CenterTitle();\r
- graph->GetXaxis()->SetNoExponent(1) ;\r
- graph->GetXaxis()->SetMoreLogLabels(1) ;\r
- graph->GetYaxis()->CenterTitle();\r
- \r
- graph->SetLineColor(color);\r
- graph->SetMarkerColor(color);\r
- graph->SetLineStyle(kDashed);\r
- graph->SetLineWidth(linewidth);\r
-\r
- return graph;\r
-\r
-}\r
-\r
-TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {\r
- //\r
- // particle = 0 ... choosen particle (setted particleMass)\r
- // particle = 1 ... Pion\r
- // particle = 2 ... Kaon\r
- // particle = 3 ... D0\r
- //\r
- Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
- \r
- Double_t particleEfficiency[400]; // with chosen particle mass\r
- Double_t kaonEfficiency[400], pionEfficiency[400], d0efficiency[400]; \r
- Double_t partEfficiency[2][400];\r
- \r
- if (particle != 0) {\r
- // resulting Pion and Kaon efficiency scaled with overall efficiency\r
- Double_t doNotDecayFactor;\r
- for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon\r
- \r
- for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
- // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.\r
- Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity\r
- if ( massloop == 1 ) { // KAON\r
- doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.\r
- kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;\r
- } else { // PION\r
- doNotDecayFactor = 1.0 ;\r
- pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ; \r
- }\r
- partEfficiency[0][j] = pionEfficiency[j];\r
- partEfficiency[1][j] = kaonEfficiency[j];\r
- } \r
- }\r
- \r
- // resulting estimate of the D0 efficiency\r
- for ( Int_t j = 0 ; j < 400 ; j++ ) {\r
- d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);\r
- }\r
- } else { \r
- for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
- particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;\r
- // NOTE: Decay factor (see kaon) should be included to be realiable\r
- }\r
- }\r
-\r
- for ( Int_t j = 0 ; j < 400 ; j++ ) { \r
- pionEfficiency[j] *= 100;\r
- kaonEfficiency[j] *= 100;\r
- d0efficiency[j] *= 100;\r
- particleEfficiency[j] *= 100;\r
- }\r
- \r
- TGraph * graph = 0;\r
- if (particle==0) {\r
- graph = new TGraph ( 400, fTransMomenta, particleEfficiency ) ; // choosen mass\r
- graph->SetLineWidth(1);\r
- } else if (particle==1) {\r
- graph = new TGraph ( 400, fTransMomenta, pionEfficiency ) ;\r
- graph->SetLineWidth(1);\r
- } else if (particle ==2) {\r
- graph = new TGraph ( 400, fTransMomenta, kaonEfficiency ) ;\r
- graph->SetLineWidth(1);\r
- } else if (particle ==3) {\r
- graph = new TGraph ( 400, fTransMomenta, d0efficiency ) ;\r
- graph->SetLineStyle(kDashed);\r
- } else \r
- return 0;\r
-\r
- graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
- graph->GetXaxis()->CenterTitle();\r
- graph->GetXaxis()->SetNoExponent(1) ;\r
- graph->GetXaxis()->SetMoreLogLabels(1) ;\r
- graph->GetYaxis()->SetTitle("Efficiency (%)") ;\r
- graph->GetYaxis()->CenterTitle();\r
- \r
- graph->SetMinimum(0.01) ; \r
- graph->SetMaximum(100) ; \r
-\r
- graph->SetLineColor(color);\r
- graph->SetMarkerColor(color);\r
- graph->SetLineWidth(linewidth);\r
-\r
- return graph;\r
-}\r
-\r
-\r
-TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {\r
- //\r
- // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)\r
- // mode 0: impact parameter (convolution of pointing and vertex resolution)\r
- // mode 1: pointing resolution\r
- // mode 2: vtx resolution \r
- \r
- \r
- TGraph *graph = new TGraph();\r
-\r
- // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?\r
- TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); // \r
- TFormula vtxResZ("vtxResZ","600/(x+6)+10"); // \r
- \r
- TGraph *trackRes = GetGraphPointingResolution(axis,1);\r
- Double_t *pt = trackRes->GetX();\r
- Double_t *trRes = trackRes->GetY();\r
- for (Int_t ip =0; ip<trackRes->GetN(); ip++) {\r
- Double_t vtxRes = 0;\r
- if (axis==0) \r
- vtxRes = vtxResRPhi.Eval(pt[ip]);\r
- else \r
- vtxRes = vtxResZ.Eval(pt[ip]);\r
- \r
- if (mode==0)\r
- graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));\r
- else if (mode ==1)\r
- graph->SetPoint(ip,pt[ip],trRes[ip]);\r
- else\r
- graph->SetPoint(ip,pt[ip],vtxRes);\r
- }\r
- \r
- graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;\r
- graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;\r
- \r
- graph->SetMinimum(1) ;\r
- graph->SetMaximum(300.1) ;\r
- graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;\r
- graph->GetXaxis()->CenterTitle();\r
- graph->GetXaxis()->SetNoExponent(1) ;\r
- graph->GetXaxis()->SetMoreLogLabels(1) ;\r
- graph->GetYaxis()->CenterTitle();\r
- \r
- graph->SetLineColor(color);\r
- graph->SetMarkerColor(color);\r
- graph->SetLineWidth(linewidth);\r
-\r
- return graph;\r
-\r
-}\r
-\r
-TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {\r
- // \r
- // returns graph according to the number\r
- //\r
- switch(number) {\r
- case 1:\r
- return GetGraphPointingResolution(0,color, linewidth); // dr\r
- case 2:\r
- return GetGraphPointingResolution(1,color, linewidth); // dz\r
- case 3:\r
- return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele\r
- case 4:\r
- return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele\r
- case 5:\r
- return GetGraphMomentumResolution(color, linewidth); // pt resolution\r
- case 10:\r
- return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle\r
- case 11:\r
- return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion\r
- case 12:\r
- return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon\r
- case 13: \r
- return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0\r
- default:\r
- printf(" Error: chosen graph number not valid\n");\r
- }\r
- return 0;\r
-\r
-}\r
-\r
-void DetectorK::MakeAliceAllNew(Bool_t flagTPC) {\r
- \r
- // All New configuration with X0 = 0.3 and resolution = 4 microns\r
- \r
- AddLayer((char*)"bpipe",2.0,0.0022); // beam pipe\r
- AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation\r
-\r
- // new ideal Pixel properties?\r
- Double_t x0 = 0.0030;\r
- Double_t resRPhi = 0.0004;\r
- Double_t resZ = 0.0004;\r
- \r
- AddLayer((char*)"ddd1", 2.2 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd2", 3.8 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd3", 6.8 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd4", 12.4 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd5", 23.5 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd6", 39.6 , x0, resRPhi, resZ); \r
- AddLayer((char*)"ddd7", 43.0 , x0, resRPhi, resZ); \r
- \r
- if (flagTPC) {\r
- AddTPC(0.1,0.1); // TPC\r
- }\r
-}\r
-\r
-void DetectorK::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) {\r
-\r
- // Numbers taken from \r
- // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks\r
- // number for misalingment: private communication with Andrea Dainese\r
-\r
- AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe\r
- AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation\r
- AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2\r
- AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2\r
-\r
-\r
- if (flagTPC) {\r
- AddTPC(0.1,0.1); // TPC\r
- }\r
- // Adding the ITS - current configuration\r
- \r
- if (AlignResiduals==0) {\r
-\r
- AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130);\r
- AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130);\r
- AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025);\r
- AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025);\r
- AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830);\r
- AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830);\r
-\r
- } else if (AlignResiduals==1) {\r
-\r
- // tracking errors ...\r
- // (Additional systematic errors due to misalignments) ... \r
- // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]\r
- // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);\r
-\r
- AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010), \r
- TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
- AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),\r
- TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
- AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), \r
- TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));\r
- AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),\r
- TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); \r
- \r
- } else if (AlignResiduals==2) {\r
- \r
- // tracking errors ... PLUS ... module misalignment\r
- \r
- // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]\r
- // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);\r
- \r
- // the ITS modules are misalignment with small gaussian smearings with\r
- // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD\r
- \r
- AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008), \r
- TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
- AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),\r
- TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));\r
- AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010), \r
- TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));\r
- AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),\r
- TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); \r
-\r
- } else {\r
- \r
- // the ITS modules are misalignment with small gaussian smearings with\r
- // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD\r
- // unknown in Z ????\r
-\r
- AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), \r
- TMath::Sqrt(0.0130*0.0130+0.000*0.000));\r
- AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),\r
- TMath::Sqrt(0.0130*0.0130+0.000*0.000));\r
- AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),\r
- TMath::Sqrt(0.0025*0.0025+0.000*0.000));\r
- AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),\r
- TMath::Sqrt(0.0025*0.0025+0.000*0.000));\r
- AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), \r
- TMath::Sqrt(0.0830*0.0830+0.000*0.000));\r
- AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),\r
- TMath::Sqrt(0.0830*0.0830+0.000*0.000)); \r
- \r
- \r
- }\r
- \r
-}\r
-\r
-\r
-void DetectorK::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {\r
- //\r
- // Produces the standard performace plots\r
- //\r
- \r
- if (!add) {\r
-\r
- TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500); \r
- c1->Divide(2,2);\r
- \r
- c1->cd(1); gPad->SetGridx(); gPad->SetGridy(); \r
- gPad->SetLogx(); \r
- TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);\r
- eff->SetTitle("Efficiencies");\r
- eff->Draw("AL");\r
- if (!onlyPionEff) {\r
- GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");\r
- GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");\r
- }\r
- c1->cd(2); gPad->SetGridx(); gPad->SetGridy(); \r
- gPad->SetLogy(); gPad->SetLogx(); \r
- GetGraphMomentumResolution(color,linewidth)->Draw("AL");\r
- \r
- c1->cd(3); gPad->SetGridx(); gPad->SetGridy(); \r
- gPad->SetLogx(); \r
- GetGraphPointingResolution(0,color,linewidth)->Draw("AL");\r
- \r
- c1->cd(4); gPad->SetGridx(); gPad->SetGridy(); \r
- gPad->SetLogx(); \r
- GetGraphPointingResolution(1,color,linewidth)->Draw("AL");\r
-\r
- } else {\r
-\r
- TVirtualPad *c1 = gPad->GetMother();\r
-\r
- c1->cd(1);\r
- GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");\r
- if (!onlyPionEff) {\r
- GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");\r
- GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");\r
- }\r
- c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");\r
- \r
- c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");\r
- \r
- c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");\r
- \r
- }\r
-\r
-}\r
-\r
-\r
-Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir) const\r
-{\r
- // Get local X of the track position estimated at the radius lab radius r. \r
- // The track curvature is accounted exactly\r
- //\r
- // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)\r
- // 0 - take the intersection closest to the current track position\r
- // >0 - go along the track (increasing fX)\r
- // <0 - go backward (decreasing fX)\r
- //\r
- // special case of R=0\r
- if (r<kAlmost0) {x=0; return kTRUE;}\r
-\r
- const double* pars = tr->GetParameter();\r
- const Double_t &fy=pars[0], &sn = pars[2];\r
- //\r
- double fx = tr->GetX();\r
- double crv = tr->GetC(bz);\r
- if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track\r
- if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis\r
- double det = (r-fx)*(r+fx);\r
- if (det<0) return kFALSE; // does not reach raduis r\r
- x = fx;\r
- if (dir==0) return kTRUE;\r
- det = TMath::Sqrt(det);\r
- if (dir>0) { // along the track direction\r
- if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle\r
- else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle\r
- }\r
- else if(dir>0) { // agains track direction\r
- if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis\r
- else if (fy>det) return kFALSE; // track is against Y axis\r
- }\r
- }\r
- else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis\r
- double det = (r-fy)*(r+fy);\r
- if (det<0) return kFALSE; // does not reach raduis r\r
- det = TMath::Sqrt(det);\r
- if (!dir) {\r
- x = fx>0 ? det : -det; // choose the solution requiring the smalest step\r
- return kTRUE;\r
- }\r
- else if (dir>0) { // along the track direction\r
- if (fx > det) return kFALSE; // current point is in on the right from the circle\r
- else if (fx <-det) x = -det; // on the left\r
- else x = det; // within the circle\r
- }\r
- else { // against the track direction\r
- if (fx <-det) return kFALSE; \r
- else if (fx > det) x = det;\r
- else x = -det;\r
- }\r
- }\r
- else { // general case of straight line\r
- double cs = TMath::Sqrt((1-sn)*(1+sn));\r
- double xsyc = fx*sn-fy*cs;\r
- double det = (r-xsyc)*(r+xsyc);\r
- if (det<0) return kFALSE; // does not reach raduis r\r
- det = TMath::Sqrt(det);\r
- double xcys = fx*cs+fy*sn;\r
- double t = -xcys;\r
- if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step\r
- else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)\r
- if (t>=-det) t += -det; // take minimal step giving t>0\r
- else return kFALSE; // both solutions have negative t\r
- }\r
- else { // go in increasing fx direction. (t+-det < 0)\r
- if (t<det) t -= det; // take minimal step giving t<0\r
- else return kFALSE; // both solutions have positive t\r
- }\r
- x = fx + cs*t;\r
- }\r
- }\r
- else { // helix\r
- // get center of the track circle\r
- double tR = 1./crv; // track radius (for the moment signed)\r
- double cs = TMath::Sqrt((1-sn)*(1+sn));\r
- double x0 = fx - sn*tR;\r
- double y0 = fy + cs*tR;\r
- double r0 = TMath::Sqrt(x0*x0+y0*y0);\r
- // printf("Xc:%+e Yc:%+e\n",x0,y0);\r
- //\r
- if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle\r
- tR = TMath::Abs(tR);\r
- double tR2r0 = tR/r0;\r
- double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);\r
- double det = (1.-g)*(1.+g);\r
- if (det<0) return kFALSE; // does not reach raduis r\r
- det = TMath::Sqrt(det);\r
- //\r
- // the intersection happens in 2 points: {x0+tR*C,y0+tR*S} \r
- // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det\r
- // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)\r
- //\r
- double tmp = 1.+g*tR2r0;\r
- x = x0*tmp; \r
- double y = y0*tmp;\r
- if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique\r
- double dfx = tR2r0*TMath::Abs(y0)*det;\r
- double dfy = tR2r0*x0*TMath::Sign(det,y0);\r
- if (dir==0) { // chose the one which corresponds to smallest step \r
- double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0\r
- if (delta<0) x += dfx;\r
- else x -= dfx;\r
- }\r
- else if (dir>0) { // along track direction: x must be > fx\r
- x -= dfx; // try the smallest step (dfx is positive)\r
- if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;\r
- }\r
- else { // backward: x must be < fx\r
- x += dfx; // try the smallest step (dfx is positive)\r
- if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;\r
- }\r
- }\r
- else { // special case: track touching the circle just in 1 point\r
- if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE; \r
- }\r
- }\r
- //\r
- return kTRUE;\r
-}\r
+#include "DetectorK.h"
+#include <TMath.h>
+#include <TMatrixD.h>
+#include <TGraph.h>
+#include <TAxis.h>
+#include <TFormula.h>
+#include <TCanvas.h>
+#include <TEllipse.h>
+#include <TText.h>
+#include <TGraphErrors.h>
+
+#include "AliExternalTrackParam.h"
+
+/***********************************************************
+
+Fast Simulation tool for Inner Tracker Systems
+
+original code of using the billoir technique was developed
+for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov
+http://rnc.lbl.gov/~jhthomas
+
+Changes by S. Rossegger -> see header file
+
+***********************************************************/
+
+
+#define RIDICULOUS 999999 // A ridiculously large resolution (cm) to flag a dead detector
+
+#define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 )
+#define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm)
+#define dNdEtaMinB 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)
+// #define dNdEtaCent 2300//15000 //1600//2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known)
+
+#define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns)
+#define AcceptanceOfTpcAndSi 1 //1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors
+#define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0)
+#define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0)
+#define EfficiencySearchFlag 2 // Define search method:
+ // -> ChiSquarePlusConfLevel = 2, ChiSquare = 1, Simple = 0.
+
+#define PionMass 0.139 // Mass of the Pion
+#define KaonMass 0.498 // Mass of the Kaon
+#define D0Mass 1.865 // Mass of the D0
+
+//TMatrixD *probKomb; // table for efficiency kombinatorics
+
+
+class CylLayerK : public TNamed {
+public:
+
+ CylLayerK(char *name) : TNamed(name,name) {}
+
+ Float_t GetRadius() const {return radius;}
+ Float_t GetRadL() const {return radL;}
+ Float_t GetPhiRes() const {return phiRes;}
+ Float_t GetZRes() const {return zRes;}
+ Float_t GetLayerEff() const {return eff;}
+
+ // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
+ Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes;
+ Float_t eff;
+ Bool_t isDead;
+
+ ClassDef(CylLayerK,1);
+};
+
+
+class ForwardLayer : public TNamed {
+public:
+ ForwardLayer(char *name) : TNamed(name,name) {}
+
+ Float_t GetZ() const {return zPos;}
+ Float_t GetXRes() const {return xRes;}
+ Float_t GetYRes() const {return yRes;}
+ Float_t GetThickness() const {return thickness;}
+ Float_t Getdensity() const {return density;}
+ Float_t GetLayerEff() const {return eff;}
+
+ // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); }
+ Float_t zPos; Float_t xRes; Float_t yRes;
+ Float_t radL;
+ Float_t thickness;
+ Float_t density;
+ Float_t eff;
+ Bool_t isDead;
+
+ ClassDef(ForwardLayer,1);
+};
+
+
+ClassImp(DetectorK)
+DetectorK::DetectorK()
+ : TNamed("test_detector","detector"),
+ fNumberOfLayers(0),
+ fNumberOfActiveLayers(0),
+ fNumberOfActiveITSLayers(0),
+ fBField(0.5),
+ fLhcUPCscale(1.0),
+ fIntegrationTime(0.02), // in ms
+ fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
+ fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
+ fParticleMass(0.140), // Standard: pion mass
+ fMaxRadiusSlowDet(10.),
+ fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
+ fAtLeastFake(1), // if at least x fakes, track is considered fake ...
+ fMaxSeedRadius(50000),
+ fptScale(10.),
+ fdNdEtaCent(2300)
+{
+ //
+ // default constructor
+ //
+ // fLayers = new TObjArray();
+
+}
+
+DetectorK::DetectorK(char *name, char *title)
+ : TNamed(name,title),
+ fNumberOfLayers(0),
+ fNumberOfActiveLayers(0),
+ fNumberOfActiveITSLayers(0),
+ fBField(0.5),
+ fLhcUPCscale(1.0),
+ fIntegrationTime(0.02), // in ms
+ fConfLevel(0.0027), // 0.27 % -> 3 sigma confidence
+ fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle
+ fParticleMass(0.140), // Standard: pion mass
+ fMaxRadiusSlowDet(10.),
+ fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers
+ fAtLeastFake(1), // if at least x fakes, track is considered fake ...
+ fMaxSeedRadius(50000),
+ fptScale(10.),
+ fdNdEtaCent(2200)
+{
+ //
+ // default constructor, that set the name and title
+ //
+ // fLayers = new TObjArray();
+}
+DetectorK::~DetectorK() { //
+ // virtual destructor
+ //
+ // delete fLayers;
+}
+
+void DetectorK::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes, Float_t eff) {
+ //
+ // Add additional layer to the list of layers (ordered by radius)
+ //
+
+ CylLayerK *newLayer = (CylLayerK*) fLayers.FindObject(name);
+
+ if (!newLayer) {
+ newLayer = new CylLayerK(name);
+ newLayer->radius = radius;
+ newLayer->radL = radL;
+ newLayer->phiRes = phiRes;
+ newLayer->zRes = zRes;
+ newLayer->eff = eff;
+
+ if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS)
+ newLayer->isDead = kTRUE;
+ else
+ newLayer->isDead = kFALSE;
+
+ if (fLayers.GetEntries()==0)
+ fLayers.Add(newLayer);
+ else {
+
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ CylLayerK *l = (CylLayerK*)fLayers.At(i);
+ if (radius<l->radius) {
+ fLayers.AddBefore(l,newLayer);
+ break;
+ }
+ if (radius>l->radius && (i+1)==fLayers.GetEntries() ) {
+ // even bigger then last one
+ fLayers.Add(newLayer);
+ }
+ }
+
+ }
+ fNumberOfLayers += 1;
+ if (!(newLayer->isDead)) {
+ fNumberOfActiveLayers += 1;
+ TString lname(newLayer->GetName());
+ if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
+ }
+
+
+ } else {
+ printf("Layer with the name %s does already exist\n",name);
+ }
+
+
+}
+
+void DetectorK::KillLayer(char *name) {
+ //
+ // Marks layer as dead. Contribution only by Material Budget
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot mark as dead\n",name);
+ else {
+ tmp->phiRes = 999999;
+ tmp->zRes = 999999;
+ if (!(tmp->isDead)) {
+ tmp->isDead = kTRUE;
+ fNumberOfActiveLayers -= 1;
+ TString lname(tmp->GetName());
+ if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
+ }
+ }
+}
+
+void DetectorK::SetRadius(char *name, Float_t radius) {
+ //
+ // Set layer radius [cm]
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+
+
+ if (!tmp) {
+ printf("Layer %s not found - cannot set radius\n",name);
+ } else {
+
+ Float_t tmpRadL = tmp->radL;
+ Float_t tmpPhiRes = tmp->phiRes;
+ Float_t tmpZRes = tmp->zRes;
+
+ RemoveLayer(name); // so that the ordering is correct
+ AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes);
+ }
+}
+
+Float_t DetectorK::GetRadius(char *name) {
+ //
+ // Return layer radius [cm]
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get radius\n",name);
+ else
+ return tmp->radius;
+
+ return 0;
+}
+
+void DetectorK::SetRadiationLength(char *name, Float_t radL) {
+ //
+ // Set layer material [cm]
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set layer material\n",name);
+ else {
+ tmp->radL = radL;
+ }
+}
+
+Float_t DetectorK::GetRadiationLength(char *name) {
+ //
+ // Return layer radius [cm]
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get layer material\n",name);
+ else
+ return tmp->radL;
+
+ return 0;
+
+}
+
+void DetectorK::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
+ //
+ // Set layer resolution in [cm]
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set resolution\n",name);
+ else {
+
+ Bool_t wasDead = tmp->isDead;
+
+ tmp->phiRes = phiRes;
+ tmp->zRes = zRes;
+ TString lname(tmp->GetName());
+
+ if (zRes==RIDICULOUS && phiRes==RIDICULOUS) {
+ tmp->isDead = kTRUE;
+ if (!wasDead) {
+ fNumberOfActiveLayers -= 1;
+ if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
+ }
+ } else {
+ tmp->isDead = kFALSE;
+ if (wasDead) {
+ fNumberOfActiveLayers += 1;
+ if (!lname.Contains("tpc")) fNumberOfActiveITSLayers += 1;
+ }
+ }
+
+
+ }
+}
+
+Float_t DetectorK::GetResolution(char *name, Int_t axis) {
+ //
+ // Return layer resolution in [cm]
+ // axis = 0: resolution in rphi
+ // axis = 1: resolution in z
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get resolution\n",name);
+ else {
+ if (axis==0) return tmp->phiRes;
+ if (axis==1) return tmp->zRes;
+ printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
+ }
+ return 0;
+}
+
+void DetectorK::SetLayerEfficiency(char *name, Float_t eff) {
+ //
+ // Set layer efficnecy (prop that his is missed within this layer)
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set layer efficiency\n",name);
+ else {
+ tmp->eff = eff;
+ }
+}
+
+Float_t DetectorK::GetLayerEfficiency(char *name) {
+ //
+ // Get layer efficnecy (prop that his is missed within this layer)
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get layer efficneicy\n",name);
+ else
+ return tmp->eff;
+
+ return 0;
+
+}
+
+void DetectorK::RemoveLayer(char *name) {
+ //
+ // Removes a layer from the list
+ //
+
+ CylLayerK *tmp = (CylLayerK*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot remove it\n",name);
+ else {
+ Bool_t wasDead = tmp->isDead;
+ fLayers.Remove(tmp);
+ fNumberOfLayers -= 1;
+ if (!wasDead) {
+ fNumberOfActiveLayers -= 1;
+ TString lname(tmp->GetName());
+ if (!lname.Contains("tpc")) fNumberOfActiveITSLayers -= 1;
+
+ }
+ }
+}
+
+
+void DetectorK::PrintLayout() {
+ //
+ // Prints the detector layout
+ //
+
+ printf("Detector %s: \"%s\"\n",GetName(),GetTitle());
+
+ if (fLayers.GetEntries()>0)
+ printf(" Name \t\t r [cm] \t X0 \t phi & z res [um] layerEff \n");
+
+ CylLayerK *tmp = 0;
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ tmp = (CylLayerK*)fLayers.At(i);
+
+ // don't print all the tpc layers
+ TString name(tmp->GetName());
+ if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
+
+ printf("%d. %s \t %03.2f \t%1.4f\t ",i,
+ tmp->GetName(), tmp->radius, tmp->radL);
+ if (tmp->phiRes==RIDICULOUS)
+ printf(" - ");
+ else
+ printf("%3.0f ",tmp->phiRes*10000);
+ if (tmp->zRes==RIDICULOUS)
+ printf(" -");
+ else
+ printf("%3.0f",tmp->zRes*10000);
+
+ if (tmp->zRes==RIDICULOUS)
+ printf("\t -\n");
+ else
+ printf("\t%0.2f\n",tmp->eff);
+
+ }
+}
+
+void DetectorK::PlotLayout(Int_t plotDead) {
+ //
+ // Plots the detector layout in Front view
+ //
+
+ Double_t x0=0, y0=0;
+
+ TGraphErrors *gr = new TGraphErrors();
+ gr->SetPoint(0,0,0);
+ CylLayerK *lastLayer = (CylLayerK*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->radius;
+ gr->SetPointError(0,maxRad,maxRad);
+ gr->Draw("APE");
+
+
+ CylLayerK *tmp = 0;
+ for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
+ tmp = (CylLayerK*)fLayers.At(i);
+
+
+ Double_t txtpos = tmp->radius;
+ if ((tmp->isDead)) txtpos*=-1; //
+ TText *txt = new TText(x0,txtpos,tmp->GetName());
+ txt->SetTextSizePixels(5); txt->SetTextAlign(21);
+ if (!tmp->isDead || plotDead) txt->Draw();
+
+ TEllipse *layEl = new TEllipse(x0,y0,tmp->radius);
+ // layEl->SetFillColor(5);
+ layEl->SetFillStyle(5001);
+ layEl->SetLineStyle(tmp->isDead+1); // dashed if not active
+ layEl->SetLineColor(4);
+ TString name(tmp->GetName());
+ if (!tmp->isDead) layEl->SetLineWidth(2);
+ if (name.Contains("tpc") ) layEl->SetLineColor(29);
+
+ if (!tmp->isDead || plotDead) layEl->Draw();
+
+ }
+
+}
+
+
+
+void DetectorK::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) {
+ //
+ // Emulates the TPC
+ //
+ // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow
+
+
+ AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage
+ AddLayer((char*)"OFC", 254.0,0.01367); // Outer Field cage
+
+ // % Radiation Lengths ... Average per TPC row (i.e. total/159 )
+ Float_t radLBoubdary = 0.0165;
+ Float_t rBoundary = 70.0; // cm
+
+ Float_t radLPerRow = 0.000036;
+
+ Float_t tpcInnerRadialPitch = 0.75 ; // cm
+ Float_t tpcMiddleRadialPitch = 1.0 ; // cm
+ Float_t tpcOuterRadialPitch = 1.5 ; // cm
+ // Float_t tpcInnerPadWidth = 0.4 ; // cm
+ // Float_t tpcMiddlePadWidth = 0.6 ; // cm
+ // Float_t tpcOuterPadWidth = 0.6 ; // cm
+ Float_t innerRows = 63 ;
+ Float_t middleRows = 64 ;
+ Float_t outerRows = 32 ;
+ Float_t tpcRows = (innerRows + middleRows + outerRows) ;
+ Float_t rowOneRadius = 85.2 ; // cm
+ Float_t row64Radius = 135.1 ; // cm
+ Float_t row128Radius = 199.2 ; // cm
+
+ // add boundary between ITS and TPC
+ AddLayer("tpc_boundary",rBoundary,radLBoubdary); // dummy errors
+
+ for ( Int_t k = 0 ; k < tpcRows ; k++ ) {
+
+ Float_t rowRadius =0;
+ if (k<innerRows)
+ rowRadius = rowOneRadius + k*tpcInnerRadialPitch ;
+ else if ( k>=innerRows && k<(innerRows+middleRows) )
+ rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ;
+ else if (k>=(innerRows+middleRows) && k<tpcRows )
+ rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ;
+
+ if ( k%skip == 0 )
+ AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean);
+ else
+ AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row
+
+
+ }
+
+}
+
+void DetectorK::RemoveTPC() {
+
+ // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
+ CylLayerK *tmp = 0;
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ tmp = (CylLayerK*)fLayers.At(i);
+ TString name(tmp->GetName());
+ if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
+ }
+ RemoveLayer((char*)"IFC");
+
+}
+
+
+Double_t DetectorK::ThetaMCS ( Double_t mass, Double_t radLength, Double_t momentum ) const
+{
+ //
+ // returns the Multiple Couloumb scattering angle (compare PDG boolet, 2010, equ. 27.14)
+ //
+
+ Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ;
+ Double_t theta = 0.0 ; // Momentum and mass in GeV
+ // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum );
+ if ( radLength > 0 ) theta = 0.0136 * TMath::Sqrt(radLength) / ( beta * momentum ) * (1+0.038*TMath::Log(radLength)) ;
+ return (theta) ;
+}
+
+
+Double_t DetectorK::ProbGoodHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
+{
+ // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm
+ // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm
+ // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius
+ Double_t sx, sy, goodHit ;
+ sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusRPhi * HitDensity(radius) ;
+ sy = 2 * TMath::Pi() * searchRadiusZ * searchRadiusZ * HitDensity(radius) ;
+ goodHit = TMath::Sqrt(1./((1+sx)*(1+sy))) ;
+ return ( goodHit ) ;
+}
+
+
+Double_t DetectorK::ProbGoodChiSqHit ( Double_t radius, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
+{
+ // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm
+ // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
+ Double_t sx, goodHit ;
+ sx = 2 * TMath::Pi() * searchRadiusRPhi * searchRadiusZ * HitDensity(radius) ;
+ goodHit = 1./(1+sx) ;
+ return ( goodHit ) ;
+}
+
+Double_t DetectorK::ProbGoodChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
+{
+ // Based on work by Ruben Shahoyen
+ // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function
+ // Plus, in addition, taking a "confidence level" and the "layer efficiency" into account
+ // Following is correct for 2 DOF
+
+ Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
+ Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
+ Double_t goodHit = leff/(2*alpha) * (1 - TMath::Exp(-alpha*c));
+ return ( goodHit ) ;
+}
+
+Double_t DetectorK::ProbNullChiSqPlusConfHit ( Double_t radius, Double_t leff, Double_t searchRadiusRPhi, Double_t searchRadiusZ )
+{
+ // Based on work by Ruben Shahoyen
+ // This is the probability to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)
+
+ Double_t c = -2 *TMath::Log(fConfLevel); // quantile at cut of confidence level
+ Double_t alpha = (1 + 2 * TMath::Pi() * HitDensity(radius) * searchRadiusRPhi * searchRadiusZ)/2;
+ Double_t nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));
+ return ( nullHit ) ;
+}
+
+Double_t DetectorK::HitDensity ( Double_t radius )
+{
+ // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor.
+ // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results
+ // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda.
+ // Note that this function assumes we are working in CM and CM**2 [not meters].
+ // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2.
+
+ // Double_t MaxRadiusSlowDet = 0.1; //? // Maximum radius for slow detectors. Fast detectors
+ // and only fast detectors reside outside this radius.
+ Double_t arealDensity = 0 ;
+
+ if ( radius > fMaxRadiusSlowDet )
+ {
+ arealDensity = OneEventHitDensity(fdNdEtaCent,radius) ; // Fast detectors see central collision density (only)
+ arealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,radius) ; // Increase density due to background
+ }
+
+ if (radius < fMaxRadiusSlowDet )
+ { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero.
+ arealDensity = OneEventHitDensity(fdNdEtaCent,radius)
+ + IntegratedHitDensity(dNdEtaMinB,radius)
+ + UpcHitDensity(radius) ;
+ arealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,radius) ;
+ // Increase density due to background
+ }
+
+ return ( arealDensity ) ;
+}
+
+
+double DetectorK::OneEventHitDensity( Double_t multiplicity, Double_t radius ) const
+{
+ // This is for one event at the vertex. No smearing.
+
+ double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
+ double tg = TMath::Tan(2*TMath::ATan(TMath::Exp(-fAvgRapidity)));
+ den = den/TMath::Sqrt(1 + 1/(tg*tg));
+
+ // double den = multiplicity / (2.*TMath::Pi()*radius*radius) ; // 2 eta ?
+ // note: surface of sphere is '4*pi*r^2'
+ // surface of cylinder is '2*pi*r* h'
+
+
+
+ return den ;
+}
+
+
+double DetectorK::IntegratedHitDensity(Double_t multiplicity, Double_t radius)
+{
+ // The integral of minBias events smeared over a gaussian vertex distribution.
+ // Based on work by Yan Lu 12/20/2006, all radii in centimeters.
+
+ Double_t zdcHz = Luminosity * 1.e-24 * CrossSectionMinB ;
+ Double_t den = zdcHz * fIntegrationTime/1000. * multiplicity * Dist(0., radius) / (2.*TMath::Pi()*radius) ;
+
+ // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex.
+ if ( den < OneEventHitDensity(multiplicity,radius) ) den = OneEventHitDensity(multiplicity,radius) ;
+
+ return den ;
+}
+
+
+double DetectorK::UpcHitDensity(Double_t radius)
+{
+ // QED electrons ...
+
+ Double_t mUPCelectrons ; ;
+ // mUPCelectrons = fLhcUPCscale * (1.23 - radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC
+ mUPCelectrons = fLhcUPCscale*5456/(radius*radius)/dNdEtaMinB; // Fit to 'Rossegger,Sadovsky'-Alice simulation
+ if ( mUPCelectrons < 0 ) mUPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R
+ mUPCelectrons *= IntegratedHitDensity(dNdEtaMinB,radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate
+ mUPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC
+
+ return mUPCelectrons ;
+}
+
+
+double DetectorK::Dist(double z, double r)
+{
+ // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution
+ // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm
+ // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters.
+ Int_t index = 1 ; // Start weight at 1 for Simpsons rule integration
+ Int_t nsteps = 301 ; // NSteps must be odd for Simpson's rule to work
+ double dist = 0.0 ;
+ double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (nsteps-1) ; //cm
+ double z0 = 0.0 ; //cm
+ for(int i=0; i<nsteps; i++){
+ if ( i == nsteps-1 ) index = 1 ;
+ z0 = -4*SigmaD + i*dz0 ;
+ dist += index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) *
+ (1/sqrt((z-z0)*(z-z0) + r*r)) ;
+ if ( index != 4 ) index = 4; else index = 2 ;
+ }
+ return dist;
+}
+
+#define PZero 0.861 // Momentum of back to back decay particles in the CM frame
+#define EPiZero 0.872 // Energy of the pion from a D0 decay at rest
+#define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest
+
+Double_t DetectorK::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) const {
+ // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi.
+
+ Double_t const1 = pt / D0Mass ;
+ Double_t const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ;
+ Double_t sum, ptPi, ptK ;
+ Double_t effp, effk ;
+
+ sum = 0.0 ;
+ for ( Int_t k = 0 ; k < 360 ; k++ ) {
+
+ Double_t theta = k * TMath::Pi() / 180. ;
+
+ ptPi = TMath::Sqrt(
+ PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
+ const1*const1*EPiZero*EPiZero -
+ 2*PZero*TMath::Cos(theta)*const2*const1*EPiZero +
+ PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
+ ) ;
+
+ ptK = TMath::Sqrt(
+ PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*const2*const2 +
+ const1*const1*EKZero*EKZero +
+ 2*PZero*TMath::Cos(theta)*const2*const1*EKZero +
+ PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta)
+ ) ;
+
+ // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays
+ Int_t pionindex = (int)((ptPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
+ Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ;
+
+ if ( pionindex >= kNptBins ) pionindex = 399 ;
+ if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ;
+ if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd
+ if ( effp < 0 ) effp = 0 ;
+
+ if ( kaonindex >= kNptBins ) kaonindex = 399 ;
+ if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ;
+ if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd
+ if ( effk < 0 ) effk = 0 ;
+
+ // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here.
+
+ sum += effp * effk ;
+
+ }
+
+ Double_t mean =sum/360;
+ return mean ;
+
+}
+
+
+
+void DetectorK::SolveViaBilloir(Int_t flagD0,Int_t print, Bool_t allPt, Double_t meanPt, char* detLayer) {
+ //
+ // Solves the current geometry with the Billoir technique
+ // ( see P. Billoir, Nucl. Instr. and Meth. 225 (1984), p. 352. )
+ // ABOVE IS OBSOLETE -> NOW, its uses the Aliroot Kalman technique
+ //
+ const float kTrackingMargin = 0.1;
+
+ static AliExternalTrackParam probTr; // track to propagate
+ probTr.SetUseLogTermMS(kTRUE);
+
+
+ Int_t nPt = kNptBins;
+ // Clean up ......
+ for (Int_t j=0; j<nPt; j++) {
+ for (Int_t i=0; i<kMaxNumberOfDetectors; i++) {
+ fDetPointRes[i][j] = RIDICULOUS;
+ fDetPointZRes[i][j] = RIDICULOUS;
+ }
+ fTransMomenta[j] =0;
+ fMomentumRes[j] =0;
+ fResolutionRPhi[j] =0;
+ }
+
+ if (!allPt) { // not the whole pt range -> allows a faster minimization at a defined 'meanpt'
+ nPt = 3;
+ }
+
+
+ // Calculate track parameters using Billoirs method of matrices
+
+ Double_t pt,tgl, pz, lambda, deltaPoverP ;
+ Double_t charge ;
+ Double_t mass[3] ;
+ Int_t printOnce = 1 ;
+
+ mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k
+
+ mass[2] = fParticleMass; // third loop
+
+ Int_t mStart =0;
+ if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode
+
+
+
+
+ // Prepare Probability Kombinations
+ Int_t nLayer = fNumberOfActiveITSLayers;
+ Int_t base = 3; // null, fake, correct
+
+ Int_t komb = (Int_t) TMath::Power(base,nLayer);
+
+ TMatrixD probLay(base,fNumberOfActiveITSLayers);
+ TMatrixD probKomb(komb,nLayer);
+ for (Int_t num=0; num<komb; num++) {
+ for (Int_t l=nLayer; l--;) {
+ Int_t pow = ((Int_t)TMath::Power(base,l+1));
+ probKomb(num,nLayer-1-l)=(num%pow)/((Int_t)TMath::Power(base,l));
+ }
+ }
+
+ TString detLayerStr(detLayer);
+ CylLayerK *theLayer = (CylLayerK*) fLayers.FindObject(detLayer);
+ if (!theLayer && detLayerStr.IsNull()!=1){
+ printf("Error: Layer with the name \"%s\" not found -> no detailed infos possible\n",detLayer);
+ return;
+ }
+ for (Int_t i=0; i<fLayers.GetEntries();i++) {
+ CylLayerK *l = (CylLayerK*) fLayers.At(i);
+ if (detLayerStr.CompareTo(l->GetName())==0) { // is the same
+ kDetLayer=i;
+ break;
+ }
+ }
+
+
+ for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) {
+
+ // PseudoRapidity OK, used as an angle
+ lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ;
+
+
+ for ( Int_t i = 0 ; i < nPt ; i++ ) { // pt loop
+
+ CylLayerK *last = (CylLayerK*) fLayers.At((fLayers.GetEntries()-1));
+
+ // Starting values based on radius of outermost layer ... log10 steps to ~20 GeV
+ Double_t bigRad = last->radius/2 ; // min. pt which the algorithm below could handle
+ // if (bigRad<61) bigRad=61; // -> min pt around 100 MeV for Bz=0.5T (don't overdo it ... ;-) )
+ fTransMomenta[i] = ( 0.3*bigRad*TMath::Abs(fBField)*1e-2 ) - 0.08 - (1./fptScale-0.1) + TMath::Power(10,2.3*i/nPt) / fptScale ;
+ if (!allPt) { // just 3 points around meanPt
+ fTransMomenta[i] = meanPt-0.001+(Double_t)(i)*0.001;
+ }
+
+ // New from here ................
+
+ // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
+ // These are the EndPoint values for y, z, a, b, and d
+ double bGauss = fBField*10; // field in kgauss
+ pt = fTransMomenta[i]; // GeV/c
+ tgl = TMath::Tan(lambda); // dip
+ charge = -1; // Assume an electron
+ pz = pt * TMath::Tan(lambda) ; // GeV/
+ enum {kY,kZ,kSnp,kTgl,kPtI}; // track parameter aliases
+ enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
+ //
+ probTr.Reset();
+ double *trPars = (double*)probTr.GetParameter();
+ double *trCov = (double*)probTr.GetCovariance();
+ trPars[kY] = 0; // start from Y = 0
+ trPars[kZ] = 0; // Z = 0
+ trPars[kSnp] = 0; // track along X axis at the vertex
+ trPars[kTgl] = TMath::Tan(lambda); // dip
+ trPars[kPtI] = charge/pt; // q/pt
+ //
+ // put tiny errors to propagate to the outer radius
+ trCov[kY2] = trCov[kZ2] = trCov[kSnp2] = trCov[kTgl2] = trCov[kPtI2] = 1e-9;
+ if (!PropagateToR(&probTr,last->radius+ kTrackingMargin,bGauss,1)) continue;
+ //if (!probTr.PropagateTo(last->radius,bGauss)) continue;
+ // reset cov.matrix
+ const double kLargeErr2Coord = 5*5;
+ const double kLargeErr2Dir = 0.7*0.7;
+ const double kLargeErr2PtI = 30.5*30.5;
+ for (int ic=15;ic--;) trCov[ic] = 0.;
+ trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
+ trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
+ trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
+ probTr.CheckCovariance();
+ //
+ // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,fTransMomenta[i],(last->radius)/100,momentum, d);
+
+ // Set Detector-Efficiency Storage area to unity
+ fEfficiency[massloop][i] = 1.0 ;
+ //
+ // Back-propagate the covariance matrix along the track.
+
+ CylLayerK *layer = 0;
+
+ // find last "active layer" - start tracking at the last active layer
+ Int_t lastActiveLayer = 0;
+ for (Int_t j=fLayers.GetEntries(); j--;) {
+ layer = (CylLayerK*)fLayers.At(j);
+ if (!(layer->isDead)) { // is alive
+ lastActiveLayer = j;
+ break;
+ }
+ }
+ // probTr.Print();
+ for (Int_t j=lastActiveLayer+1; j--;) { // Layer loop
+
+ layer = (CylLayerK*)fLayers.At(j);
+
+ if (layer->radius>fMaxSeedRadius) continue; // no seeding beyond this radius
+
+ TString name(layer->GetName());
+ Bool_t isVertex = name.Contains("vertex");
+ //
+ if (!PropagateToR(&probTr,layer->radius,bGauss,-1)) exit(1);
+ // if (!probTr.PropagateTo(last->radius,bGauss)) exit(1); //
+ // rotate to frame with X axis normal to the surface
+ if (!isVertex) {
+ double pos[3];
+ probTr.GetXYZ(pos); // lab position
+ double phi = TMath::ATan2(pos[1],pos[0]);
+ if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
+ if (!probTr.Rotate(phi)) {
+ printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f (pt=%+.3f)\n",
+ phi,layer->radius,pos[0],pos[1],pos[2],pt);
+
+ probTr.Print();
+ exit(1);
+ }
+ }
+ // save resolutions at this layer
+ fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
+ fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
+ //printf(">> L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
+ // End save
+ //
+ if (isVertex) continue;
+ //
+ // create fake measurement with the errors assigned to the layer
+ // account for the measurement there
+ double meas[2] = {probTr.GetY(),probTr.GetZ()};
+ double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
+ //
+
+ if (!probTr.Update(meas,measErr2)) {
+ printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
+ meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
+ probTr.Print();
+ exit(1);
+ }
+
+ // correct for materials of this layer
+ // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
+ if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
+ printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
+ probTr.Print();
+ exit(1);
+ }
+ //
+ }
+
+ // Pattern recognition is done .... save values like vertex resolution etc.
+
+ // Convert the Convariance matrix parameters into physical quantities
+ // The results are propogated to the previous point but *do not* include the measurement at that point.
+ // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/probTr.Get1P(); // Absolute magnitude so ignore charge
+ deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/TMath::Abs(probTr.GetSigned1Pt());
+ fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent
+ fResolutionRPhi[i] = TMath::Sqrt( probTr.GetSigmaY2() ) * 1.e4; // result in microns
+ fResolutionZ[i] = TMath::Sqrt( probTr.GetSigmaZ2() ) * 1.e4; // result in microns
+ // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
+
+
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+ printf("Number of active layers: %d\n",fNumberOfActiveLayers) ;
+ if (fAtLeastCorr != -1) printf("Number of combinatorics for probabilities: %d\n",komb);
+ printf("Mass of tracked particle: %f (at pt=%5.0lf MeV)\n",fParticleMass,fTransMomenta[i]*1000);
+ printf("Name Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ;
+ // printOnce =0;
+ }
+
+ // print out and efficiency calculation
+ Int_t iLayActive=0;
+ for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
+
+ layer = (CylLayerK*)fLayers.At(j);
+
+ // Convert to Meters, Tesla, and GeV
+ Float_t radius = layer->radius /100;
+ Float_t phiRes = layer->phiRes /100;
+ Float_t zRes = layer->zRes /100;
+ Float_t radLength = layer->radL;
+ Float_t leff = layer->eff; // basic layer efficiency
+ Bool_t isDead = layer->isDead;
+
+
+ if ( (!isDead && radLength >0) ) {
+
+ Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
+ phiRes * phiRes ) * 100. ; // work in cm
+ Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
+ zRes * zRes ) * 100. ; // work in cm
+
+
+ Double_t layerEfficiency = 0;
+ if ( EfficiencySearchFlag == 0 )
+ layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
+ else if ( EfficiencySearchFlag == 1 )
+ layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
+ else if ( EfficiencySearchFlag == 2 )
+ layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
+
+ TString name(layer->GetName());
+ if (!name.Contains("tpc")) {
+ probLay(2,iLayActive)= layerEfficiency ; // Pcorr
+ probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
+ probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
+ iLayActive++;
+ }
+ if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
+
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+
+
+ printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
+ layer->GetName(), radius*100, radLength,
+ fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
+ phiRes*1.e6, zRes*1.e6,
+ HitDensity(radius*100)) ;
+ if (!name.Contains("tpc"))
+ printf("%10.3f\n", layerEfficiency);
+ else
+ printf(" - \n");
+ }
+
+ if (!name.Contains("tpc") && !name.Contains("trd")) fEfficiency[massloop][i] *= layerEfficiency;
+
+
+ }
+
+ if (fAtLeastCorr != -1) {
+ // Calculate probabilities from Kombinatorics tree ...
+ Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
+ fEfficiency[massloop][i] = probs[0]; // efficiency
+ fFake[massloop][i] = probs[1]; // fake
+ }
+
+ /*
+ // vertex print
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1 && radius==0) {
+ printf("%s:\t ----- ----- %10.0f %11.0f \n", layer->GetName(),fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6);
+ }
+ */
+ }
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+ if (fNumberOfActiveLayers >=1500) printOnce = 0 ;
+ printf("\n") ;
+ }
+
+
+
+
+ if (fNumberOfActiveLayers <1500 ) {
+
+
+
+ // BACKWORD TRACKING +++++++++++++++++
+ // number of layers is quite low ... efficiency calculation was probably nonsense
+ // Tracking outward (backword) to get reliable efficiencies from "smoothed estimates"
+
+ // For below, see paper, NIM A262 (1987) p.444, eqs.12.
+ // Equivalently, one can simply combine the forward and backward estimates. Assuming
+ // pf,Cf and pb,Cb as extrapolated position estimates and errors from fwd and bwd passes one can
+ // use a weighted estimate Cw = (Cf^-1 + Cb^-1)^-1, pw = Cw (pf Cf^-1 + pb Cb^-1).
+ // Surely, for the most extreme point, where one error matrices is infinite, this does not change anything.
+
+ Bool_t doLikeAliRoot = 0; // don't do the "combined info" but do like in Aliroot
+
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+ printf("- Numbers of active layer is low (%d):\n -> \"outward\" fitting done as well to get reliable eff.estimates\n",
+ fNumberOfActiveLayers);
+ }
+
+ // RESET Covariance Matrix ( to 10 x the estimate -> as it is done in AliExternalTrackParam)
+ // mIstar.UnitMatrix(); // start with unity
+ if (doLikeAliRoot) {
+ probTr.ResetCovariance(10);
+ } else {
+ // cannot do complete reset, set to very large errors
+ for (int ic=15;ic--;) trCov[ic] = 0.;
+ trCov[kY2] = trCov[kZ2] = kLargeErr2Coord;
+ trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
+ trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
+ probTr.CheckCovariance();
+ // cout<<pt<<": "<<kLargeErr2Coord<<" "<<kLargeErr2Dir<<" "<<kLargeErr2PtI*trPars[kPtI]*trPars[kPtI]<<endl;
+ }
+ // Clean up and storing of "forward estimates"
+ Double_t detPointResForw[kMaxNumberOfDetectors][kNptBins], detPointZResForw[kMaxNumberOfDetectors][kNptBins] ;
+ for (Int_t k=0; k<kMaxNumberOfDetectors; k++) {
+ for (Int_t l=0; l<nPt; l++) {
+ detPointResForw[k][l] = fDetPointRes[k][l];
+ if (!doLikeAliRoot) fDetPointRes[k][l] = RIDICULOUS;
+ detPointZResForw[k][l] = fDetPointZRes[k][l];
+ if (!doLikeAliRoot) fDetPointZRes[k][l] = RIDICULOUS;
+ }
+ }
+
+ // find first "active layer" - start tracking at the first active layer
+ Int_t firstActiveLayer = 0;
+ for (Int_t j=0; j<fLayers.GetEntries(); j++) {
+ layer = (CylLayerK*)fLayers.At(j);
+ if (!(layer->isDead)) { // is alive
+ firstActiveLayer = j;
+ break;
+ }
+ }
+ probTr.Rotate(0);
+ for (Int_t j=firstActiveLayer; j<(fLayers.GetEntries()); j++) { // Layer loop
+
+ layer = (CylLayerK*)fLayers.At(j);
+ // CylLayerK *nextlayer = (CylLayerK*)fLayers.At(j+1);
+
+ TString name(layer->GetName());
+ Bool_t isVertex = name.Contains("vertex");
+ if (!PropagateToR(&probTr, layer->radius,bGauss,1)) exit(1);
+ //if (!probTr.PropagateTo(last->radius,bGauss)) exit(1);
+ if (!isVertex) {
+ // rotate to frame with X axis normal to the surface
+ double pos[3];
+ probTr.GetXYZ(pos); // lab position
+ double phi = TMath::ATan2(pos[1],pos[0]);
+ if ( TMath::Abs(TMath::Abs(phi)-TMath::Pi()/2)<1e-3) phi = 0;//TMath::Sign(TMath::Pi()/2 - 1e-3,phi);
+ if (!probTr.Rotate(phi)) {
+ printf("Failed to rotate to the frame (phi:%+.3f)of layer at %.2f at XYZ: %+.3f %+.3f %+.3f (pt=%+.3f)\n",
+ phi,layer->radius,pos[0],pos[1],pos[2],pt);
+
+ probTr.Print();
+ exit(1);
+ }
+ }
+ //
+ fDetPointRes [j][i] = TMath::Sqrt( probTr.GetSigmaY2() )/100 ; // result in meters
+ fDetPointZRes[j][i] = TMath::Sqrt( probTr.GetSigmaZ2() )/100 ; // result in meters
+ //
+ //printf("<< L%d r:%e sy: %e sz: %e\n",j,layer->radius,fDetPointRes[j][i],fDetPointZRes[j][i]);
+ // create fake measurement with the errors assigned to the layer
+ // account for the measurement there
+ if (isVertex) continue;
+ double meas[2] = {probTr.GetY(),probTr.GetZ()};
+ double measErr2[3] = {layer->phiRes*layer->phiRes,0,layer->zRes*layer->zRes};
+ //
+ if (!probTr.Update(meas,measErr2)) {
+ printf("Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}\n",
+ meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]);
+ probTr.Print();
+ exit(1);
+ }
+ // correct for materials of this layer
+ // note: if apart from MS we want also e.loss correction, the density*length should be provided as 2nd param
+ if (!probTr.CorrectForMeanMaterial(layer->radL, 0, mass[massloop] , kTRUE)) {
+ printf("Failed to apply material correction, X/X0=%.4f\n",layer->radL);
+ probTr.Print();
+ exit(1);
+ }
+ }
+
+ // values below NOT REALIABLE -> they do not point to the vertex but outwards !!!!!!!
+ // ++++++++++++++
+ // also update the values for the track position ??????
+ /*
+ // Pattern recognition is done .... save values like vertex resolution etc.
+
+ // Invert the Matrix to recover the convariance matrix
+ mIstar.Invert() ;
+ // Convert the Convariance matrix parameters into physical quantities
+ // The results are propogated to the previous point but *do not* include the measurement at that point.
+ deltaPoverP = TMath::Sqrt( mIstar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore charge
+ fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ) ; // results in percent
+ fResolutionRPhi[i] = TMath::Sqrt( mIstar(0,0) ) * 1.e6 ; // result in microns
+ fResolutionZ[i] = TMath::Sqrt( mIstar(1,1) ) * 1.e6 ; // result in microns
+ // equivalent[i] = TMath::Sqrt(fResolutionRPhi[i]*fResolutionZ[i]) ; // Equivalent circular radius
+ */
+
+ // deltaPoverP = TMath::Sqrt(probTr.GetSigma1Pt2())/TMath::Abs(probTr.GetSigned1Pt());
+ // fMomentumRes[i] = 100.* TMath::Abs( deltaPoverP ); // results in percent
+
+ // Weighted combination of the forward and backward estimates
+ if (!doLikeAliRoot) {
+ for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) {
+ fDetPointRes[j][i] = detPointResForw[j][i]*fDetPointRes[j][i]/TMath::Sqrt((detPointResForw[j][i]*detPointResForw[j][i]) + (fDetPointRes[j][i]*fDetPointRes[j][i]));
+ fDetPointZRes[j][i] = detPointZResForw[j][i]*fDetPointZRes[j][i]/TMath::Sqrt((detPointZResForw[j][i]*detPointZResForw[j][i]) + (fDetPointZRes[j][i]*fDetPointZRes[j][i]));
+ }
+ }
+ // Set Detector-Efficiency Storage area to unity
+ fEfficiency[massloop][i] = 1.0 ;
+
+ // print out and efficiency calculation
+ iLayActive=0;
+ for (Int_t j=(fLayers.GetEntries()-1); j>=0; j--) { // Layer loop
+
+ layer = (CylLayerK*)fLayers.At(j);
+
+ // Convert to Meters, Tesla, and GeV
+ Float_t radius = layer->radius /100;
+ Float_t phiRes = layer->phiRes /100;
+ Float_t zRes = layer->zRes /100;
+ Float_t radLength = layer->radL;
+ Float_t leff = layer->eff;
+ Bool_t isDead = layer->isDead;
+
+
+ Double_t layerEfficiency = 0;
+ if ( (!isDead && radLength >0) ) {
+ Double_t rphiError = TMath::Sqrt( fDetPointRes[j][i] * fDetPointRes [j][i] +
+ phiRes * phiRes ) * 100. ; // work in cm
+ Double_t zError = TMath::Sqrt( fDetPointZRes[j][i] * fDetPointZRes[j][i] +
+ zRes * zRes ) * 100. ; // work in cm
+ if ( EfficiencySearchFlag == 0 )
+ layerEfficiency = ProbGoodHit( radius*100, rphiError , zError ) ;
+ else if ( EfficiencySearchFlag == 1 )
+ layerEfficiency = ProbGoodChiSqHit( radius*100, rphiError , zError ) ;
+ else if ( EfficiencySearchFlag == 2 )
+ layerEfficiency = ProbGoodChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ;
+
+ TString name(layer->GetName());
+ if (!name.Contains("tpc")) {
+ probLay(2,iLayActive)= layerEfficiency ; // Pcorr
+ probLay(0,iLayActive)= ProbNullChiSqPlusConfHit( radius*100,leff, rphiError , zError ) ; // Pnull
+ probLay(1,iLayActive)= 1 - probLay(2,iLayActive) - probLay(0,iLayActive); // Pfake
+ iLayActive++;
+ }
+ if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue;
+
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+
+
+ printf("%s:\t%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f ",
+ layer->GetName(), radius*100, radLength,
+ fDetPointRes[j][i]*1.e6, fDetPointZRes[j][i]*1.e6,
+ phiRes*1.e6, zRes*1.e6,
+ HitDensity(radius*100)) ;
+ if (!name.Contains("tpc"))
+ printf("%10.3f\n", layerEfficiency);
+ else
+ printf(" - \n");
+ }
+
+ if (massloop==2 && j==kDetLayer) { // copy layer specific performances
+ fEfficProlongLay[i] = layerEfficiency;
+ }
+
+ if (!name.Contains("tpc") && !name.Contains("trd")) fEfficiency[massloop][i] *= layerEfficiency;
+
+
+
+ }
+ if (fAtLeastCorr != -1) {
+ // Calculate probabilities from Kombinatorics tree ...
+ Double_t *probs = PrepareEffFakeKombinations(&probKomb, &probLay);
+ fEfficiency[massloop][i] = probs[0]; // efficiency
+ fFake[massloop][i] = probs[1]; // fake
+ }
+ }
+ if (print == 1 && fTransMomenta[i] >= meanPt && massloop == 2 && printOnce == 1) {
+ printOnce = 0 ;
+ printf("\n") ;
+ }
+ }
+
+ if (massloop==2) { // copy layer specific performances
+ fResolutionRPhiLay[i] = fDetPointRes[kDetLayer][i];
+ fResolutionZLay[i] = fDetPointZRes[kDetLayer][i];
+ }
+
+ } // pt loop
+
+
+
+ } // mass loop
+
+ probTr.SetUseLogTermMS(kFALSE); // Reset of MS term usage to avoid problems since its static
+
+
+
+}
+
+
+TGraph * DetectorK::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
+ //
+ // returns the momentum resolution
+ //
+
+ TGraph *graph = new TGraph(kNptBins, fTransMomenta, fMomentumRes);
+ graph->SetTitle("Momentum Resolution .vs. Pt" ) ;
+ // graph->GetXaxis()->SetRangeUser(0.,5.0) ;
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->SetTitle("Momentum Resolution (%)") ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetMaximum(20);
+ graph->SetMinimum(0.1);
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+
+}
+
+TGraph * DetectorK::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) {
+
+ // Returns the pointing resolution
+ // axis = 0 ... rphi pointing resolution
+ // axis = 1 ... z pointing resolution
+ //
+
+ TGraph * graph = 0;
+
+ if (axis==0) {
+ graph = new TGraph ( kNptBins, fTransMomenta, fResolutionRPhi ) ;
+ graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
+ } else {
+ graph = new TGraph ( kNptBins, fTransMomenta, fResolutionZ ) ;
+ graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
+ }
+
+ graph->SetMinimum(1) ;
+ graph->SetMaximum(1000.1) ;
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetLineWidth(linewidth);
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+
+ return graph;
+
+}
+
+TGraph * DetectorK::GetGraphLayerInfo(Int_t plot, Int_t color, Int_t linewidth) {
+
+ // Returns the pointing resolution
+ // plot = 0 ... rphi pointing resolution
+ // plot = 1 ... z pointing resolution
+ // plot = 2 ... prolongation efficiency (outwards)
+ //
+
+
+ Double_t fDet[kNptBins];
+ for ( Int_t i = 0 ; i < kNptBins ; i++ ) { // pt loop
+ if (plot==0)
+ fDet[i] = fResolutionRPhiLay[i]*1e6; // in microns
+ else if (plot==1)
+ fDet[i] = fResolutionZLay[i]*1e6; // in microns
+ else
+ fDet[i] = fEfficProlongLay[i]*100; // in percent
+ }
+
+ CylLayerK *l = (CylLayerK*) fLayers.At(kDetLayer);
+ TGraph * graph = 0;
+ graph = new TGraph ( kNptBins, fTransMomenta, fDet ) ;
+ if (plot==0) {
+ graph->SetTitle(Form("R-#phi Pointing Resolution onto layer \"%s\"",(char*)l->GetName()) );
+ graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
+ } else if (plot==1){
+ graph->SetTitle(Form("Z Pointing Resolution onto layer \"%s\"",(char*)l->GetName()) ) ;
+ graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum)") ;
+ } else {
+ graph->SetTitle(Form("Prolongation efficiency onto layer \"%s\"",(char*)l->GetName()) ) ;
+ graph->GetYaxis()->SetTitle("Prolongation efficiency (%)") ;
+ graph->SetMinimum(0);
+ graph->SetMaximum(100);
+ }
+
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetLineWidth(linewidth);
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+
+ return graph;
+
+}
+
+
+
+TGraph * DetectorK::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) {
+ //
+ // returns the Pointing resolution (accoring to Telescope equation)
+ // axis =0 ... in rphi
+ // axis =1 ... in z
+ //
+
+ Double_t resolution[kNptBins];
+
+ Double_t layerResolution[2];
+ Double_t layerRadius[2];
+ Double_t layerThickness[2];
+
+ Int_t count =0; // search two first active layers
+ printf("Telescope equation for layers: ");
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ CylLayerK *l = (CylLayerK*)fLayers.At(i);
+ if (!l->isDead && l->radius>0) {
+ layerRadius[count] = l->radius;
+ layerThickness[count] = l->radL;
+ if (axis==0) {
+ layerResolution[count] = l->phiRes;
+ } else {
+ layerResolution[count] = l->zRes;
+ }
+ printf("%s, ",l->GetName());
+ count++;
+ }
+ if (count>=2) break;
+ }
+ printf("\n");
+
+ Double_t pt, momentum, thickness,aMCS ;
+ Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
+
+ for ( Int_t i = 0 ; i < kNptBins ; i++ ) {
+ // Reference data as if first two layers were acting all alone
+ pt = fTransMomenta[i] ;
+ momentum = pt / TMath::Cos(lambda) ; // Total momentum
+ resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1]
+ + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ;
+ resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ;
+ thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ;
+ aMCS = ThetaMCS(fParticleMass, thickness, momentum) ;
+ resolution[i] += layerRadius[0]*layerRadius[0]*aMCS*aMCS ;
+ resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns
+ }
+
+
+
+ TGraph* graph = new TGraph ( kNptBins, fTransMomenta, resolution ) ;
+
+ if (axis==0) {
+ graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("RPhi Pointing Resolution (#mum) ") ;
+ } else {
+ graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("Z Pointing Resolution (#mum) ") ;
+ }
+ graph->SetMinimum(1) ;
+ graph->SetMaximum(300.1) ;
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineStyle(kDashed);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+
+}
+
+TGraph * DetectorK::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) {
+ //
+ // particle = 0 ... choosen particle (setted particleMass)
+ // particle = 1 ... Pion
+ // particle = 2 ... Kaon
+ // particle = 3 ... D0
+ //
+ Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
+
+ Double_t particleEfficiency[kNptBins]; // with chosen particle mass
+ Double_t kaonEfficiency[kNptBins], pionEfficiency[kNptBins], d0efficiency[kNptBins];
+ Double_t partEfficiency[2][400];
+
+ if (particle != 0) {
+ // resulting Pion and Kaon efficiency scaled with overall efficiency
+ Double_t doNotDecayFactor;
+ for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
+
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
+ Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
+ if ( massloop == 1 ) { // KAON
+ doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
+ kaonEfficiency[j] = fEfficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
+ } else { // PION
+ doNotDecayFactor = 1.0 ;
+ pionEfficiency[j] = fEfficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ;
+ }
+ partEfficiency[0][j] = pionEfficiency[j];
+ partEfficiency[1][j] = kaonEfficiency[j];
+ }
+ }
+
+ // resulting estimate of the D0 efficiency
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
+ }
+ } else {
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
+ // NOTE: Decay factor (see kaon) should be included to be realiable
+ }
+ }
+
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ pionEfficiency[j] *= 100;
+ kaonEfficiency[j] *= 100;
+ d0efficiency[j] *= 100;
+ particleEfficiency[j] *= 100;
+ }
+
+ TGraph * graph = 0;
+ if (particle==0) {
+ graph = new TGraph ( kNptBins, fTransMomenta, particleEfficiency ) ; // choosen mass
+ graph->SetLineWidth(1);
+ } else if (particle==1) {
+ graph = new TGraph ( kNptBins, fTransMomenta, pionEfficiency ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==2) {
+ graph = new TGraph ( kNptBins, fTransMomenta, kaonEfficiency ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==3) {
+ graph = new TGraph ( kNptBins, fTransMomenta, d0efficiency ) ;
+ graph->SetLineStyle(kDashed);
+ } else
+ return 0;
+
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->SetTitle("Efficiency (%)") ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetMinimum(0.01) ;
+ graph->SetMaximum(100) ;
+
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+}
+
+TGraph * DetectorK::GetGraphRecoFakes(Int_t particle,Int_t color, Int_t linewidth) {
+ //
+ // particle = 0 ... choosen particle (setted particleMass)
+ // particle = 1 ... Pion
+ // particle = 2 ... Kaon
+ //
+
+ Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
+
+ Double_t particleFake[kNptBins]; // with chosen particle mass
+ Double_t kaonFake[kNptBins], pionFake[kNptBins];
+ Double_t partFake[2][kNptBins];
+
+ if (particle != 0) {
+ // resulting Pion and Kaon efficiency scaled with overall efficiency
+ Double_t doNotDecayFactor;
+ for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon
+
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm.
+ Double_t momentum = fTransMomenta[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity
+ if ( massloop == 1 ) { // KAON
+ doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm.
+ kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;
+ } else { // PION
+ pionFake[j] = fFake[0][j] ;
+ }
+ partFake[0][j] = pionFake[j];
+ partFake[1][j] = kaonFake[j];
+ }
+ }
+
+ } else {
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ particleFake[j] = fFake[2][j];
+ // NOTE: Decay factor (see kaon) should be included to be realiable
+ }
+ }
+
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ pionFake[j] *= 100;
+ kaonFake[j] *= 100;
+ particleFake[j] *= 100;
+ }
+
+ TGraph * graph = 0;
+ if (particle==0) {
+ graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
+ graph->SetLineWidth(1);
+ } else if (particle==1) {
+ graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==2) {
+ graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
+ graph->SetLineWidth(1);
+ }
+
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->SetTitle("Fake (%)") ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetMinimum(0.01) ;
+ graph->SetMaximum(100) ;
+
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+}
+TGraph * DetectorK::GetGraphRecoPurity(Int_t particle,Int_t color, Int_t linewidth) {
+ //
+ // particle = 0 ... choosen particle (setted particleMass)
+ // particle = 1 ... Pion
+ // particle = 2 ... Kaon
+ //
+
+ // Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity));
+
+ Double_t particleFake[kNptBins]; // with chosen particle mass
+ Double_t kaonFake[kNptBins], pionFake[kNptBins];
+ // Double_t partFake[2][kNptBins];
+
+ if (particle != 0) {
+ cout <<" not implemented"<<endl;
+
+ } else {
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ particleFake[j] = fFake[2][j];
+ // NOTE: Decay factor (see kaon) should be included to be realiable
+ }
+ }
+
+ // Get Purity
+ for ( Int_t j = 0 ; j < kNptBins ; j++ ) {
+ pionFake[j] = (1-pionFake[j])*100;
+ kaonFake[j] = (1-kaonFake[j])*100;
+ particleFake[j] = (1-particleFake[j])*100;
+ }
+
+ TGraph * graph = 0;
+ if (particle==0) {
+ graph = new TGraph ( kNptBins, fTransMomenta, particleFake ) ; // choosen mass
+ graph->SetLineWidth(1);
+ } else if (particle==1) {
+ graph = new TGraph ( kNptBins, fTransMomenta, pionFake ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==2) {
+ graph = new TGraph ( kNptBins, fTransMomenta, kaonFake ) ;
+ graph->SetLineWidth(1);
+ }
+
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->SetTitle("Purity (%)") ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetMinimum(0.01) ;
+ graph->SetMaximum(100) ;
+
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+}
+
+
+TGraph* DetectorK::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) {
+ //
+ // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution)
+ // mode 0: impact parameter (convolution of pointing and vertex resolution)
+ // mode 1: pointing resolution
+ // mode 2: vtx resolution
+
+
+ TGraph *graph = new TGraph();
+
+ // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ?
+ TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); //
+ TFormula vtxResZ("vtxResZ","600/(x+6)+10"); //
+
+ TGraph *trackRes = GetGraphPointingResolution(axis,1);
+ Double_t *pt = trackRes->GetX();
+ Double_t *trRes = trackRes->GetY();
+ for (Int_t ip =0; ip<trackRes->GetN(); ip++) {
+ Double_t vtxRes = 0;
+ if (axis==0)
+ vtxRes = vtxResRPhi.Eval(pt[ip]);
+ else
+ vtxRes = vtxResZ.Eval(pt[ip]);
+
+ if (mode==0)
+ graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip]));
+ else if (mode ==1)
+ graph->SetPoint(ip,pt[ip],trRes[ip]);
+ else
+ graph->SetPoint(ip,pt[ip],vtxRes);
+ }
+
+ graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("d_{0} r#phi resolution (#mum)") ;
+
+ graph->SetMinimum(1) ;
+ graph->SetMaximum(300.1) ;
+ graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ;
+ graph->GetXaxis()->CenterTitle();
+ graph->GetXaxis()->SetNoExponent(1) ;
+ graph->GetXaxis()->SetMoreLogLabels(1) ;
+ graph->GetYaxis()->CenterTitle();
+
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+ graph->SetLineWidth(linewidth);
+
+ return graph;
+
+}
+
+TGraph* DetectorK::GetGraph(Int_t number, Int_t color, Int_t linewidth) {
+ //
+ // returns graph according to the number
+ //
+ switch(number) {
+ case 1:
+ return GetGraphPointingResolution(0,color, linewidth); // dr
+ case 2:
+ return GetGraphPointingResolution(1,color, linewidth); // dz
+ case 3:
+ return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele
+ case 4:
+ return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele
+ case 5:
+ return GetGraphMomentumResolution(color, linewidth); // pt resolution
+ case 10:
+ return GetGraphRecoEfficiency(0, color, linewidth); // tracked particle
+ case 11:
+ return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion
+ case 12:
+ return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon
+ case 13:
+ return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0
+ case 15:
+ return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle
+ case 16:
+ return GetGraphRecoFakes(1, color, linewidth); // Fake pion
+ case 17:
+ return GetGraphRecoFakes(2, color, linewidth); // Fake kaon
+ default:
+ printf(" Error: chosen graph number not valid\n");
+ }
+ return 0;
+
+}
+
+void DetectorK::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon) {
+
+ // All New configuration with X0 = 0.3 and resolution = 4 microns
+
+ AddLayer((char*)"bpipe",2.0,0.0022); // beam pipe
+ AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
+
+ // new ideal Pixel properties?
+ Double_t x0 = 0.0050;
+ Double_t resRPhi = 0.0006;
+ Double_t resZ = 0.0006;
+
+ if (flagMon) {
+ x0 = 0.0030;
+ resRPhi = 0.0004;
+ resZ = 0.0004;
+ }
+
+ AddLayer((char*)"ddd1", 2.2 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd2", 2.8 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd3", 3.6 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd4", 20.0 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd5", 22.0 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd6", 41.0 , x0, resRPhi, resZ);
+ AddLayer((char*)"ddd7", 43.0 , x0, resRPhi, resZ);
+
+ if (flagTPC) {
+ AddTPC(0.1,0.1); // TPC
+ }
+}
+
+void DetectorK::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) {
+
+ // Numbers taken from
+ // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks
+ // number for misalingment: private communication with Andrea Dainese
+
+ AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe
+ AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation
+ AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2
+ AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2
+
+
+ if (flagTPC) {
+ AddTPC(0.1,0.1); // TPC
+ }
+ // Adding the ITS - current configuration
+
+ if (AlignResiduals==0) {
+
+ AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0130);
+ AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0130);
+ AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025);
+ AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025);
+ AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830);
+ AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830);
+
+ } else if (AlignResiduals==1) {
+
+ // tracking errors ...
+ // (Additional systematic errors due to misalignments) ...
+ // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
+ // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
+
+ AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
+ TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
+ AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
+ TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
+
+ } else if (AlignResiduals==2) {
+
+ // tracking errors ... PLUS ... module misalignment
+
+ // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm]
+ // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000);
+
+ // the ITS modules are misalignment with small gaussian smearings with
+ // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
+
+ AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
+ AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
+
+ } else {
+
+ // the ITS modules are misalignment with small gaussian smearings with
+ // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD
+ // unknown in Z ????
+
+ AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.000*0.000));
+ AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.000*0.000));
+ AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.000*0.000));
+ AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.000*0.000));
+ AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.000*0.000));
+ AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.000*0.000));
+
+
+ }
+
+}
+
+
+void DetectorK::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
+ //
+ // Produces the standard performace plots
+ //
+ TGraph *eff,*momRes,*pointRes;
+ if (!add) {
+
+ TCanvas *c1 = new TCanvas("c1","c1");//,100,100,500,500);
+ c1->Divide(2,2);
+
+ c1->cd(1); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogx();
+ eff = GetGraphRecoEfficiency(1,color,linewidth);
+ eff->SetName(Form("grEff%d",1));
+ eff->SetTitle("Efficiencies");
+ eff->Draw("AL");
+ eff->SetMaximum(110);
+ if (!onlyPionEff) {
+ eff = GetGraphRecoEfficiency(2,color,linewidth);
+ eff->SetName(Form("grEff%d",2));
+ eff->Draw("L");
+ eff = GetGraphRecoEfficiency(3,color,linewidth);
+ eff->SetName(Form("grEff%d",3));
+ eff->Draw("L");
+ }
+ c1->cd(2); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogy(); gPad->SetLogx();
+ momRes = GetGraphMomentumResolution(color,linewidth);
+ momRes->SetName(Form("grMomRes%d",1));
+ momRes->Draw("AL");
+
+ c1->cd(3); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogx();
+ pointRes = GetGraphPointingResolution(0,color,linewidth);
+ pointRes->SetName(Form("pointRRes%d",0));
+ pointRes->Draw("AL");
+ //
+ c1->cd(4); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogx();
+ pointRes = GetGraphPointingResolution(1,color,linewidth);
+ pointRes->SetName(Form("pointZRes%d",0));
+ pointRes->Draw("AL");
+
+ } else {
+
+ TVirtualPad *c1 = gPad->GetMother();
+
+ c1->cd(1);
+ eff = GetGraphRecoEfficiency(1,color,linewidth);
+ eff->SetName(Form("grEff%dadd",1));
+ eff->Draw("L");
+ if (!onlyPionEff) {
+ eff = GetGraphRecoEfficiency(2,color,linewidth);
+ eff->SetName(Form("grEff%dadd",2));
+ eff->Draw("L");
+ eff = GetGraphRecoEfficiency(3,color,linewidth);
+ eff->SetName(Form("grEff%dadd",3));
+ eff->Draw("L");
+ }
+ c1->cd(2);
+ momRes = GetGraphMomentumResolution(color,linewidth);
+ momRes->SetName(Form("grMomRes%dadd",1));
+ momRes->Draw("L");
+
+ c1->cd(3);
+ pointRes = GetGraphPointingResolution(0,color,linewidth);
+ pointRes->SetName(Form("pointRRes%dadd",0));
+ pointRes->Draw("L");
+
+ c1->cd(4);
+ pointRes = GetGraphPointingResolution(1,color,linewidth);
+ pointRes->SetName(Form("pointZRes%dadd",0));
+ pointRes->Draw("L");
+
+ }
+
+}
+
+
+Bool_t DetectorK::GetXatLabR(AliExternalTrackParam* tr,Double_t r,Double_t &x, Double_t bz, Int_t dir)
+{
+ // Get local X of the track position estimated at the radius lab radius r.
+ // The track curvature is accounted exactly
+ //
+ // The flag "dir" can be used to remove the ambiguity of which intersection to take (out of 2 possible)
+ // 0 - take the intersection closest to the current track position
+ // >0 - go along the track (increasing fX)
+ // <0 - go backward (decreasing fX)
+ //
+ // special case of R=0
+ if (r<kAlmost0) {x=0; return kTRUE;}
+
+ const double* pars = tr->GetParameter();
+ const Double_t &fy=pars[0], &sn = pars[2];
+ //
+ double fx = tr->GetX();
+ double crv = tr->GetC(bz);
+ if (TMath::Abs(crv)<=kAlmost0) { // this is a straight track
+ if (TMath::Abs(sn)>=kAlmost1) { // || to Y axis
+ double det = (r-fx)*(r+fx);
+ if (det<0) return kFALSE; // does not reach raduis r
+ x = fx;
+ if (dir==0) return kTRUE;
+ det = TMath::Sqrt(det);
+ if (dir>0) { // along the track direction
+ if (sn>0) {if (fy>det) return kFALSE;} // track is along Y axis and above the circle
+ else {if (fy<-det) return kFALSE;} // track is against Y axis amd belo the circle
+ }
+ else if(dir>0) { // agains track direction
+ if (sn>0) {if (fy<-det) return kFALSE;} // track is along Y axis
+ else if (fy>det) return kFALSE; // track is against Y axis
+ }
+ }
+ else if (TMath::Abs(sn)<=kAlmost0) { // || to X axis
+ double det = (r-fy)*(r+fy);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ if (!dir) {
+ x = fx>0 ? det : -det; // choose the solution requiring the smalest step
+ return kTRUE;
+ }
+ else if (dir>0) { // along the track direction
+ if (fx > det) return kFALSE; // current point is in on the right from the circle
+ else if (fx <-det) x = -det; // on the left
+ else x = det; // within the circle
+ }
+ else { // against the track direction
+ if (fx <-det) return kFALSE;
+ else if (fx > det) x = det;
+ else x = -det;
+ }
+ }
+ else { // general case of straight line
+ double cs = TMath::Sqrt((1-sn)*(1+sn));
+ double xsyc = fx*sn-fy*cs;
+ double det = (r-xsyc)*(r+xsyc);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ double xcys = fx*cs+fy*sn;
+ double t = -xcys;
+ if (dir==0) t += t>0 ? -det:det; // chose the solution requiring the smalest step
+ else if (dir>0) { // go in increasing fX direction. ( t+-det > 0)
+ if (t>=-det) t += -det; // take minimal step giving t>0
+ else return kFALSE; // both solutions have negative t
+ }
+ else { // go in increasing fx direction. (t+-det < 0)
+ if (t<det) t -= det; // take minimal step giving t<0
+ else return kFALSE; // both solutions have positive t
+ }
+ x = fx + cs*t;
+ }
+ }
+ else { // helix
+ // get center of the track circle
+ double tR = 1./crv; // track radius (for the moment signed)
+ double cs = TMath::Sqrt((1-sn)*(1+sn));
+ double x0 = fx - sn*tR;
+ double y0 = fy + cs*tR;
+ double r0 = TMath::Sqrt(x0*x0+y0*y0);
+ // printf("Xc:%+e Yc:%+e tR:%e r0:%e\n",x0,y0,tR,r0);
+ //
+ if (r0<=kAlmost0) return kFALSE; // the track is concentric to circle
+ tR = TMath::Abs(tR);
+ double tR2r0 = tR/r0;
+ double g = 0.5*(r*r/(r0*tR) - tR2r0 - 1./tR2r0);
+ double det = (1.-g)*(1.+g);
+ if (det<0) return kFALSE; // does not reach raduis r
+ det = TMath::Sqrt(det);
+ //
+ // the intersection happens in 2 points: {x0+tR*C,y0+tR*S}
+ // with C=f*c0+-|s0|*det and S=f*s0-+c0 sign(s0)*det
+ // where s0 and c0 make direction for the circle center (=x0/r0 and y0/r0)
+ //
+ double tmp = 1.+g*tR2r0;
+ x = x0*tmp;
+ double y = y0*tmp;
+ if (TMath::Abs(y0)>kAlmost0) { // when y0==0 the x,y is unique
+ double dfx = tR2r0*TMath::Abs(y0)*det;
+ double dfy = tR2r0*x0*TMath::Sign(det,y0);
+ if (dir==0) { // chose the one which corresponds to smallest step
+ double delta = (x-fx)*dfx-(y-fy)*dfy; // the choice of + in C will lead to smaller step if delta<0
+ if (delta<0) x += dfx;
+ else x -= dfx;
+ }
+ else if (dir>0) { // along track direction: x must be > fx
+ x -= dfx; // try the smallest step (dfx is positive)
+ if (x<fx && (x+=dfx+dfx)<fx) return kFALSE;
+ }
+ else { // backward: x must be < fx
+ x += dfx; // try the smallest step (dfx is positive)
+ if (x>fx && (x-=dfx+dfx)>fx) return kFALSE;
+ }
+ }
+ else { // special case: track touching the circle just in 1 point
+ if ( (dir>0&&x<fx) || (dir<0&&x>fx) ) return kFALSE;
+ }
+ }
+ //
+ return kTRUE;
+}
+
+
+
+Double_t* DetectorK::PrepareEffFakeKombinations(TMatrixD *probKomb, TMatrixD *probLay) {
+
+ if (!probLay) {
+ printf("Error: Layer tracking efficiencies not set \n");
+ return 0;
+ }
+
+ TMatrixD &tProbKomb = *probKomb;
+ TMatrixD &tProbLay = *probLay;
+
+
+ // Int_t base = tProbLay.GetNcols(); // 3? null, fake, correct
+ Int_t nLayer = tProbKomb.GetNcols(); // nlayer? - number of ITS layers
+ Int_t komb = tProbKomb.GetNrows(); // 3^nlayer? - number of kombinations
+
+ // Fill probabilities
+
+ Double_t probEff =0;
+ Double_t probFake =0;
+ for (Int_t num=0; num<komb; num++) {
+ Int_t flCorr=0, flFake=0, flNull=0;
+ for (Int_t l=0; l<nLayer; l++) {
+ if (tProbKomb(num,l)==0)
+ flNull++;
+ else if (tProbKomb(num,l)==1)
+ flFake++;
+ else if (tProbKomb(num,l)==2)
+ flCorr++;
+ else
+ printf("Error: unexpected values in combinatorics table\n");
+ }
+
+ Int_t fkAtLeastCorr = fAtLeastCorr;
+ if (fAtLeastCorr == -1) fkAtLeastCorr = nLayer; // all hits are "correct"
+
+ if (flCorr>=fkAtLeastCorr && flFake==0) { // at least correct but zero fake
+ Double_t probEffLayer = 1;
+ for (Int_t l=0; l<nLayer; l++) {
+ probEffLayer *= tProbLay((Int_t)tProbKomb(num,l),l);
+ // cout<<a(num,l)<<" ";
+ }
+ // cout<<endl;
+ probEff+=probEffLayer;
+ }
+
+ if (flFake>=fAtLeastFake) {
+ Double_t probFakeLayer = 1;
+ for (Int_t l=0; l<nLayer; l++) {
+ probFakeLayer *= tProbLay((Int_t)tProbKomb(num,l),l);
+ // cout<<a(num,l)<<" ";
+ }
+ // cout<<endl;
+ probFake+=probFakeLayer;
+ }
+
+ }
+ Double_t *probs = new Double_t[2];
+ probs[0] = probEff; probs[1] = probFake;
+ return probs;
+
+}
+
+//____________________________________
+Bool_t DetectorK::PropagateToR(AliExternalTrackParam* trc, double r, double b, int dir)
+{
+ // go to radius R
+ //
+ double xR = 0;
+ double rr = r*r;
+ int iter = 0;
+ const double kTiny = 1e-6;
+ while(1) {
+ if (!GetXatLabR(trc, r ,xR, b, dir)) {
+ printf("Track with pt=%f cannot reach radius %f\n",trc->Pt(),r);
+ trc->Print();
+ return kFALSE;
+ }
+ if (!trc->PropagateTo(xR, b)) {printf("Failed to propagate to X=%f for R=%f\n",xR,r); trc->Print(); return kFALSE;}
+ double rcurr2 = xR*xR + trc->GetY()*trc->GetY();
+ if (TMath::Abs(rcurr2-rr)<kTiny || rr<kTiny) return kTRUE;
+ //
+ // two radii correspond to this X...
+ double pos[3]; trc->GetXYZ(pos);
+ double phi = TMath::ATan2(pos[1],pos[0]); //TMath::ASin( trc->GetSnp() );
+ if (!trc->Rotate(phi)) {printf("Failed to rotate to %f to propagate to R=%f\n",phi,r); trc->Print(); return kFALSE;}
+ if (++iter>8) {printf("Failed to propagate to R=%f after %d steps\n",r,iter); trc->Print(); return kFALSE;}
+ }
+ return kTRUE;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff