-#include "KMCDetector.h"\r
-#include <TMath.h>\r
-#include <TH1F.h>\r
-#include <TProfile.h>\r
-#include <TF1.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 <TRandom3.h>\r
-#include <TGraphErrors.h>\r
-#include <TStopwatch.h>\r
-\r
-/***********************************************************\r
-\r
-Fast Simulation tool for Inner Tracker Systems\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 1//950//660//950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700)\r
-// #define dNdEtaCent 2300//15000 //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
-//TMatrixD *probKomb; // table for efficiency kombinatorics\r
-\r
-ClassImp(KMCProbe)\r
-\r
-Int_t KMCProbe::fgNITSLayers = 0;\r
-Double_t KMCProbe::fgMissingHitPenalty = 2.;\r
-//__________________________________________________________________________\r
-KMCProbe& KMCProbe::operator=(const KMCProbe& src) \r
-{\r
- if (this!=&src) {\r
- AliExternalTrackParam::operator=(src);\r
- fMass = src.fMass;\r
- fChi2 = src.fChi2;\r
- fHits = src.fHits;\r
- fFakes = src.fFakes;\r
- fNHits = src.fNHits;\r
- fNHitsITS = src.fNHitsITS;\r
- fNHitsITSFake = src.fNHitsITSFake;\r
- fInnLrCheck = src.fInnLrCheck;\r
- for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];\r
- }\r
- return *this;\r
-}\r
-\r
-//__________________________________________________________________________\r
-KMCProbe::KMCProbe(KMCProbe& src) \r
- : AliExternalTrackParam(src),\r
- fMass(src.fMass),\r
- fChi2(src.fChi2),\r
- fHits(src.fHits),\r
- fFakes(src.fFakes),\r
- fNHits(src.fNHits),\r
- fNHitsITS(src.fNHitsITS),\r
- fNHitsITSFake(src.fNHitsITSFake),\r
- fInnLrCheck(src.fInnLrCheck)\r
-{\r
- for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];\r
-}\r
-\r
-//__________________________________________________________________________\r
-void KMCProbe::ResetCovMat()\r
-{\r
- // reset errors\r
- double *trCov = (double*)GetCovariance();\r
- double *trPars = (double*)GetParameter();\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
-}\r
-\r
-//__________________________________________________________________________\r
-void KMCProbe::Print(Option_t* option) const\r
-{\r
- printf("M=%.3f Chi2=%7.2f (Norm:%6.2f|%d) Hits: Total:%d ITS:%d ITSFakes:%d | Y:%+8.4f Z: %+8.4f |", \r
- fMass,fChi2,GetNormChi2(kTRUE),fInnLrCheck,fNHits,fNHitsITS,fNHitsITSFake, GetY(),GetZ());\r
- for (int i=0;i<fgNITSLayers;i++) {\r
- if (!IsHit(i)) printf(".");\r
- else printf("%c",IsHitFake(i) ? '-':'+');\r
- }\r
- printf("|%s\n",IsKilled() ? " KILLED":"");\r
- TString opt = option;\r
- if (!opt.IsNull()) AliExternalTrackParam::Print(option);\r
-}\r
-\r
-//__________________________________________________________________________\r
-Int_t KMCProbe::Compare(const TObject* obj) const\r
-{\r
- // compare to tracks\r
- const KMCProbe* trc = (KMCProbe*) obj;\r
- if (trc->IsKilled()) {\r
- if (IsKilled()) return 0;\r
- return -1;\r
- }\r
- else if (IsKilled()) return 1;\r
- double chi2a = GetNormChi2(kTRUE);\r
- double chi2b = trc->GetNormChi2(kTRUE);\r
- if (chi2a<chi2b) return -1;\r
- if (chi2a>chi2b) return 1;\r
- return 0;\r
-}\r
-\r
-//__________________________________________________________________________\r
-Bool_t KMCProbe::GetXatLabR(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 = GetParameter();\r
- const Double_t &fy=pars[0], &sn = pars[2];\r
- //\r
- double fx = GetX();\r
- double crv = 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) {\r
- AliDebug(2,Form("r0 = %f",r0));\r
- return kFALSE;\r
- } // 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) {\r
- AliDebug(2,Form("g=%f tR=%f r0=%f\n",g,tR, r0));\r
- return kFALSE;\r
- } // 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
-\r
-//____________________________________\r
-Bool_t KMCProbe::PropagateToR(double r, double b, int dir) \r
-{\r
- // go to radius R\r
- //\r
- double xR = 0;\r
- double rr = r*r;\r
- int iter = 0;\r
- const double kTiny = 1e-4;\r
- while(1) {\r
- if (!GetXatLabR(r ,xR, b, dir)) {\r
- // printf("Track with pt=%f cannot reach radius %f\n",Pt(),r);\r
- // Print("l");\r
- return kFALSE;\r
- }\r
- \r
- if (!PropagateTo(xR, b)) {\r
- if (AliLog::GetGlobalDebugLevel()>2) {\r
- printf("Failed to propagate to X=%f for R=%f\n",xR,r); \r
- Print("l"); \r
- }\r
- return kFALSE;\r
- }\r
- double rcurr2 = xR*xR + GetY()*GetY();\r
- if (TMath::Abs(rcurr2-rr)<kTiny || rr<kAlmost0) return kTRUE;\r
- //\r
- // two radii correspond to this X...\r
- double pos[3]; GetXYZ(pos);\r
- double phi = TMath::ATan2(pos[1],pos[0]);\r
- if (!Rotate(phi)) {\r
- if (AliLog::GetGlobalDebugLevel()>2) {\r
- printf("Failed to rotate to %f to propagate to R=%f\n",phi,r); \r
- Print("l"); \r
- }\r
- return kFALSE;\r
- }\r
- if (++iter>8) {\r
- if (AliLog::GetGlobalDebugLevel()>2) {\r
- printf("Failed to propagate to R=%f after %d steps\n",r,iter); \r
- Print("l"); \r
- }\r
- return kFALSE;\r
- }\r
- } \r
- return kTRUE;\r
-}\r
-\r
-\r
-//__________________________________________________________________________\r
-Bool_t KMCProbe::CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward)\r
-{\r
- // printf("before at r=%.1f p=%.4f\n",lr->fR, P());\r
- if (AliExternalTrackParam::CorrectForMeanMaterial(lr->fx2X0, inward ? lr->fXRho : -lr->fXRho, GetMass() , kTRUE)) {\r
- // printf("after at r=%.1f p=%.4f\n",lr->fR, P());\r
- return kTRUE;\r
- }\r
- AliDebug(2,Form("Failed to apply material correction, X/X0=%.4f", lr->fx2X0));\r
- if (AliLog::GetGlobalDebugLevel()>1) Print();\r
- return kFALSE;\r
-}\r
-\r
-/////////////////////////////////////////////////////////////////////////////\r
-ClassImp(KMCCluster)\r
-\r
-//_________________________________________________________________________\r
-KMCCluster::KMCCluster(KMCCluster &src) \r
-: TObject(src),\r
- fY(src.fY),fZ(src.fZ),fX(src.fX),fPhi(src.fPhi)\r
-{}\r
-\r
-//__________________________________________________________________________\r
-KMCCluster& KMCCluster::operator=(const KMCCluster& src) \r
-{\r
- if (this!=&src) {\r
- TObject::operator=(src);\r
- fY = src.fY;\r
- fZ = src.fZ;\r
- fX = src.fX;\r
- fPhi = src.fPhi;\r
- }\r
- return *this;\r
-}\r
-\r
-//_________________________________________________________________________\r
-void KMCCluster::Print(Option_t *) const \r
-{\r
- printf(" Local YZ = (%3.4lf,%3.4lf) | X=%3.4lf phi: %+.3f %s\n",fY,fZ,fX,fPhi,IsKilled()?"Killed":""); \r
-}\r
-\r
-/////////////////////////////////////////////////////////////////////////////\r
-ClassImp(KMCLayer)\r
-\r
-Double_t KMCLayer::fgDefEff = 1.0;\r
-//__________________________________________________________________________\r
-KMCLayer::KMCLayer(char *name) : \r
- TNamed(name,name),fR(0),fx2X0(0),fPhiRes(0),fZRes(0),fEff(0),fIsDead(kFALSE),fType(-1),fActiveID(-1),fSig2EstD(999),fSig2EstZ(999),\r
- fClCorr(),fClMC(),fClBg("KMCCluster",5), fTrCorr(), fTrMC("KMCProbe",5)\r
-{\r
- Reset();\r
-}\r
-\r
-//__________________________________________________________________________\r
-void KMCLayer::Reset() \r
-{\r
- fTrCorr.Reset();\r
- fClCorr.Reset();\r
- ResetMC();\r
- fSig2EstD = fSig2EstZ = 999;\r
- //\r
-}\r
-\r
-//__________________________________________________________________________\r
-KMCProbe* KMCLayer::AddMCTrack(KMCProbe* src) \r
-{\r
- int ntr = GetNMCTracks(); \r
- KMCProbe* prb = 0;\r
- if (src) prb = new(fTrMC[ntr]) KMCProbe(*src);\r
- else prb = new(fTrMC[ntr]) KMCProbe();\r
- if (!IsDead()) prb->ResetHit(GetActiveID());\r
- return prb;\r
-}\r
-\r
-//__________________________________________________________________________\r
-void KMCLayer::Print(Option_t *opt) const\r
-{\r
- printf("Lr%3d(A%3d) %10s R=%5.1f X2X0=%.3f XRho=%.3f SigY=%.4f SigZ=%.4f Eff:%4.2f\n",GetUniqueID(),fActiveID,GetName(), fR, fx2X0,fXRho,fPhiRes,fZRes,fEff);\r
- TString opts = opt; opts.ToLower();\r
- if (opts.Contains("c")) {\r
- printf("Clusters: MC: %+7.4f:%+7.4f Ideal: %+7.4f:%+7.4f NBgCl: %3d NTrMC: %4d\n",fClMC.fY,fClMC.fZ, fClCorr.fY,fClCorr.fZ, GetNBgClusters(),GetNMCTracks());\r
- }\r
-}\r
-\r
-/////////////////////////////////////////////////////////////////////////////\r
-Double_t KMCDetector::fgVtxConstraint[2]={-1,-1};\r
-\r
-ClassImp(KMCDetector)\r
-KMCDetector::KMCDetector() :\r
-TNamed("test_detector","detector"),\r
- fNLayers(0),\r
- fNActiveLayers(0),\r
- fNActiveITSLayers(0),\r
- fLastActiveLayer(-1),\r
- fLastActiveITSLayer(-1),\r
- fLastActiveLayerTracked(-1),\r
- fBFieldG(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
- fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers\r
- fAtLeastFake(1), // if at least x fakes, track is considered fake ...\r
- fdNdEtaCent(2300),\r
- fMaxChi2Cl(25.),\r
- fMaxNormChi2NDF(5.),\r
- fMinITSHits(4),\r
-//\r
- fMaxChi2ClSQ(4.), // precalulated internally\r
- fMaxSeedToPropagate(300),\r
-//\r
- fUseBackground(kFALSE),\r
- fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),\r
- fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),\r
- fDensFactorEta(1.),\r
-//\r
- fUpdCalls(0),\r
- fHMCLrResidRPhi(0),\r
- fHMCLrResidZ(0),\r
- fHMCLrChi2(0),\r
- fPattITS(0)\r
-{\r
- //\r
- // default constructor\r
- //\r
- // fLayers = new TObjArray();\r
- RequireMaxChi2Cl(fMaxChi2Cl); // just to precalulate default square\r
-}\r
-\r
-KMCDetector::KMCDetector(char *name, char *title)\r
- : TNamed(name,title),\r
- fNLayers(0),\r
- fNActiveLayers(0),\r
- fNActiveITSLayers(0),\r
- fLastActiveLayer(-1),\r
- fLastActiveITSLayer(-1), \r
- fLastActiveLayerTracked(-1),\r
- fBFieldG(5.0),\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
- fAtLeastCorr(-1), // if -1, then correct hit on all ITS layers\r
- fAtLeastFake(1), // if at least x fakes, track is considered fake ...\r
- fdNdEtaCent(2300),\r
- fMaxChi2Cl(9.),\r
- fMaxNormChi2NDF(5.),\r
- fMinITSHits(4),\r
- //\r
- fMaxChi2ClSQ(3.), // precalulated internally\r
- fMaxSeedToPropagate(50),\r
- //\r
- fUseBackground(kFALSE),\r
- fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),\r
- fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),\r
- fDensFactorEta(1.),\r
- //\r
- fUpdCalls(0),\r
- fHMCLrResidRPhi(0),\r
- fHMCLrResidZ(0),\r
- fHMCLrChi2(0),\r
- fPattITS(0)\r
-{\r
- //\r
- // default constructor, that set the name and title\r
- //\r
- // fLayers = new TObjArray();\r
-}\r
-KMCDetector::~KMCDetector() { // \r
- // virtual destructor\r
- //\r
- // delete fLayers;\r
-}\r
-\r
-void KMCDetector::InitMCWatchHistos()\r
-{\r
- // init utility histos used for MC tuning\r
- enum {kNBinRes=1000};\r
- const double kMaxResidRPhi=1.0,kMaxResidZ=1.0,kMaxChi2=50;\r
- int nlr = fNActiveITSLayers;\r
- TString nam = "mc_residrhi";\r
- TString tit = "MC $Delta Cl-Tr R#phi";\r
- fHMCLrResidRPhi = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidRPhi,nlr,-0.5,nlr-0.5);\r
- fHMCLrResidRPhi->GetXaxis()->SetTitle("cl-tr #Delta r#phi");\r
- fHMCLrResidRPhi->Sumw2();\r
- //\r
- nam = "mc_residz";\r
- tit = "MC $Delta Cl-Tr Z";\r
- fHMCLrResidZ = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidZ,nlr,-0.5,nlr-0.5);\r
- fHMCLrResidZ->GetXaxis()->SetTitle("cl-tr #Delta z");\r
- fHMCLrResidZ->Sumw2();\r
- //\r
- nam = "mc_chi2";\r
- tit = "MC #chi^{2} Cl-Tr Z";\r
- fHMCLrChi2 = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxChi2,nlr,-0.5,nlr-0.5);\r
- fHMCLrChi2->GetXaxis()->SetTitle("cl-tr #chi^{2}");\r
- fHMCLrChi2->Sumw2();\r
- //\r
- SetBit(kUtilHisto);\r
-}\r
-\r
-\r
-void KMCDetector::AddLayer(char *name, Float_t radius, Float_t x2X0, Float_t xrho, 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
- KMCLayer *newLayer = (KMCLayer*) fLayers.FindObject(name);\r
-\r
- if (!newLayer) {\r
- newLayer = new KMCLayer(name);\r
- newLayer->fR = radius;\r
- newLayer->fx2X0 = x2X0;\r
- newLayer->fXRho = xrho;\r
- newLayer->fPhiRes = phiRes;\r
- newLayer->fZRes = zRes;\r
- eff = TMath::Min(1.f,eff);\r
- newLayer->fEff = eff <0 ? KMCLayer::GetDefEff() : eff;\r
- newLayer->fActiveID = -2;\r
- TString lname = name;\r
- newLayer->fType = KMCLayer::kTypeNA;\r
- if (lname.Contains("tpc")) newLayer->fType = KMCLayer::kTPC;\r
- else if (lname.Contains("its")) newLayer->fType = KMCLayer::kITS;\r
- if (lname.Contains("vertex")) newLayer->SetBit(KMCLayer::kBitVertex);\r
- //\r
- if (newLayer->fType==KMCLayer::kTypeNA) printf("Attention: the layer %s has undefined type\n",name);\r
- //\r
- newLayer->fIsDead = (newLayer->fPhiRes==RIDICULOUS && newLayer->fZRes==RIDICULOUS);\r
- //\r
- if (fLayers.GetEntries()==0) \r
- fLayers.Add(newLayer);\r
- else {\r
- //\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- KMCLayer *l = (KMCLayer*)fLayers.At(i);\r
- if (radius<l->fR) { fLayers.AddBefore(l,newLayer); break; }\r
- if (radius>l->fR && (i+1)==fLayers.GetEntries() ) fLayers.Add(newLayer); // even bigger then last one\r
- }\r
- //\r
- }\r
- //\r
- ClassifyLayers();\r
- //\r
- } else {\r
- printf("Layer with the name %s does already exist\n",name);\r
- }\r
-}\r
-\r
-//____________________________________________________________\r
-void KMCDetector::ClassifyLayers()\r
-{\r
- // assign active Id's, etc\r
- fLastActiveLayer = -1;\r
- fLastActiveITSLayer = -1;\r
- fNActiveLayers = 0;\r
- fNActiveITSLayers = 0;\r
- //\r
- int nl = GetNLayers();\r
- for (int il=0;il<nl;il++) {\r
- KMCLayer* lr = GetLayer(il);\r
- lr->SetUniqueID(il);\r
- if (!lr->IsDead()) {\r
- fLastActiveLayer = il; \r
- lr->fActiveID = fNActiveLayers++;\r
- if (lr->IsITS()) {\r
- fLastActiveITSLayer = il;\r
- fNActiveITSLayers++;\r
- }\r
- }\r
- }\r
- //\r
- KMCProbe::SetNITSLayers(fNActiveITSLayers);\r
-}\r
-\r
-//____________________________________________________________\r
-void KMCDetector::KillLayer(char *name) {\r
- //\r
- // Marks layer as dead. Contribution only by Material Budget\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot mark as dead\n",name);\r
- else {\r
- tmp->fPhiRes = 999999;\r
- tmp->fZRes = 999999;\r
- tmp->fIsDead = kTRUE;\r
- ClassifyLayers();\r
- }\r
-}\r
-\r
-void KMCDetector::SetRadius(char *name, Float_t radius) {\r
- //\r
- // Set layer radius [cm]\r
- //\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- //\r
- if (!tmp) {\r
- printf("Layer %s not found - cannot set radius\n",name);\r
- } else {\r
- \r
- Float_t tmpRadL = tmp->fx2X0;\r
- Float_t tmpPhiRes = tmp->fPhiRes;\r
- Float_t tmpZRes = tmp->fZRes;\r
- Float_t tmpXRho = tmp->fXRho;\r
- RemoveLayer(name); // so that the ordering is correct\r
- AddLayer(name,radius,tmpRadL,tmpXRho,tmpPhiRes,tmpZRes);\r
- }\r
-}\r
-\r
-Float_t KMCDetector::GetRadius(char *name) {\r
- //\r
- // Return layer radius [cm]\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get radius\n",name);\r
- else \r
- return tmp->fR;\r
-\r
- return 0;\r
-}\r
-\r
-void KMCDetector::SetRadiationLength(char *name, Float_t x2X0) {\r
- //\r
- // Set layer material [cm]\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set layer material\n",name);\r
- else {\r
- tmp->fx2X0 = x2X0;\r
- }\r
-}\r
-\r
-Float_t KMCDetector::GetRadiationLength(char *name) {\r
- //\r
- // Return layer radius [cm]\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get layer material\n",name);\r
- else \r
- return tmp->fx2X0;\r
- \r
- return 0;\r
- \r
-}\r
-\r
-void KMCDetector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {\r
- //\r
- // Set layer resolution in [cm]\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set resolution\n",name);\r
- else {\r
- tmp->fPhiRes = phiRes;\r
- tmp->fZRes = zRes;\r
- tmp->fIsDead = (zRes==RIDICULOUS && phiRes==RIDICULOUS);\r
- ClassifyLayers();\r
- }\r
-}\r
-\r
-Float_t KMCDetector::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
- KMCLayer *tmp = GetLayer(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get resolution\n",name);\r
- else {\r
- if (axis==0) return tmp->fPhiRes;\r
- if (axis==1) return tmp->fZRes;\r
- printf("error: axis must be either 0 or 1 (rphi or z axis)\n");\r
- }\r
- return 0;\r
-}\r
-\r
-void KMCDetector::SetLayerEfficiency(char *name, Float_t eff) {\r
- //\r
- // Set layer efficnecy (prop that his is missed within this layer) \r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot set layer efficiency\n",name);\r
- else {\r
- tmp->fEff = eff;\r
- }\r
-}\r
-\r
-Float_t KMCDetector::GetLayerEfficiency(char *name) {\r
- //\r
- // Get layer efficnecy (prop that his is missed within this layer) \r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot get layer efficneicy\n",name);\r
- else \r
- return tmp->fEff;\r
- \r
- return 0;\r
- \r
-}\r
-\r
-void KMCDetector::RemoveLayer(char *name) {\r
- //\r
- // Removes a layer from the list\r
- //\r
-\r
- KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);\r
- if (!tmp) \r
- printf("Layer %s not found - cannot remove it\n",name);\r
- else {\r
- fLayers.Remove(tmp);\r
- ClassifyLayers();\r
- }\r
-}\r
-\r
-\r
-void KMCDetector::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
- KMCLayer *tmp = 0;\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- tmp = (KMCLayer*)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->fR, tmp->fx2X0);\r
- if (tmp->fPhiRes==RIDICULOUS) \r
- printf(" - ");\r
- else\r
- printf("%3.0f ",tmp->fPhiRes*10000);\r
- if (tmp->fZRes==RIDICULOUS) \r
- printf(" -\n");\r
- else\r
- printf("%3.0f\n",tmp->fZRes*10000);\r
- }\r
-}\r
-\r
-void KMCDetector::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
- KMCLayer *lastLayer = (KMCLayer*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->fR;\r
- gr->SetPointError(0,maxRad,maxRad);\r
- gr->Draw("APE");\r
- \r
-\r
- KMCLayer *tmp = 0;\r
- for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {\r
- tmp = (KMCLayer*)fLayers.At(i);\r
- \r
-\r
- Double_t txtpos = tmp->fR;\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->fR);\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 KMCDetector::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*)"IFCtpc", 77.8,0.01367, 0); // Inner Field cage (RS: set correct xrho for eloss)\r
- \r
- // % Radiation Lengths ... Average per TPC row (i.e. total/159 )\r
- Float_t x2X0PerRow = 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,x2X0PerRow,0, phiResMean,zResMean); \r
- else \r
- AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0); // non "active" row\r
- \r
- \r
- }\r
- \r
-}\r
-\r
-void KMCDetector::RemoveTPC() {\r
-\r
- // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)\r
- KMCLayer *tmp = 0;\r
- for (Int_t i = 0; i<fLayers.GetEntries(); i++) {\r
- tmp = (KMCLayer*)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 KMCDetector::ThetaMCS ( Double_t mass, Double_t x2X0, 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 ( x2X0 > 0 ) theta = 0.0136 * TMath::Sqrt(x2X0) / ( beta * momentum ) * (1+0.038*TMath::Log(x2X0)) ;\r
- return (theta) ;\r
-}\r
-\r
-\r
-Double_t KMCDetector::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 KMCDetector::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 KMCDetector::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 c = -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*c));\r
- return ( goodHit ) ; \r
-}\r
-\r
-Double_t KMCDetector::ProbNullChiSqPlusConfHit ( 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 to not have any match to the track (see also :ProbGoodChiSqPlusConfHit:)\r
-\r
- Double_t c = -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 nullHit = (1-leff+fConfLevel*leff)*TMath::Exp(-c*(alpha-1./2));\r
- return ( nullHit ) ; \r
-}\r
-\r
-Double_t KMCDetector::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
- if (radius<0.01) return 0;\r
-\r
- if ( radius >= fMaxRadiusSlowDet ) \r
- {\r
- arealDensity = OneEventHitDensity(fdNdEtaCent,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(fdNdEtaCent,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 KMCDetector::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) * fDensFactorEta ; // 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 KMCDetector::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
- double dens1 = OneEventHitDensity(multiplicity,radius);\r
- if ( den < dens1 ) den = dens1;\r
-\r
- return den ;\r
-} \r
-\r
-\r
-double KMCDetector::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 KMCDetector::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 KMCDetector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) 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(fBFieldG)) ; \r
- Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ; \r
- \r
- if ( pionindex >= 20 ) 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 >= 20 ) 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
-KMCProbe* KMCDetector::PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi, double x,double y, double z)\r
-{\r
- // Prepare trackable Kalman track at the farthest position\r
- //\r
- // Set track parameters\r
- // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis\r
- fProbe.Reset();\r
- fProbe.SetMass(mass);\r
- KMCProbe* probe = new KMCProbe(fProbe);\r
- double *trPars = (double*)probe->GetParameter();\r
- double *trCov = (double*)probe->GetCovariance();\r
- double xyz[3] = {x,y,z};\r
- probe->Global2LocalPosition(xyz,phi);\r
- probe->Set(xyz[0],phi,trPars,trCov);\r
- trPars[KMCProbe::kY] = xyz[1];\r
- trPars[KMCProbe::kZ] = xyz[2];\r
- trPars[KMCProbe::kSnp] = 0; // track along X axis at the vertex\r
- trPars[KMCProbe::kTgl] = TMath::Tan(lambda); // dip\r
- trPars[KMCProbe::kPtI] = charge/pt; // q/pt \r
- //\r
- // put tiny errors to propagate to the outer-most radius\r
- trCov[KMCProbe::kY2] = trCov[KMCProbe::kZ2] = trCov[KMCProbe::kSnp2] = trCov[KMCProbe::kTgl2] = trCov[KMCProbe::kPtI2] = 1e-20;\r
- fProbe = *probe; // store original track\r
- //\r
- // propagate to last layer\r
- fLastActiveLayerTracked = 0;\r
- for (Int_t j=0; j<=fLastActiveLayer; j++) {\r
- KMCLayer* lr = GetLayer(j);\r
- lr->Reset();\r
- //\r
- if (!PropagateToLayer(probe,lr,1)) break;\r
- if (!probe->CorrectForMeanMaterial(lr, kFALSE)) break;\r
- //\r
- lr->fClCorr.Set(probe->GetY(),probe->GetZ(), probe->GetX(), probe->GetAlpha());\r
- if (!lr->IsDead()) fLastActiveLayerTracked = j;\r
- }\r
- probe->ResetCovMat();// reset cov.matrix\r
- printf("Last active layer trracked: %d (out of %d)\n",fLastActiveLayerTracked,fLastActiveLayer);\r
- //\r
- return probe;\r
-}\r
-\r
-\r
-\r
-TGraph * KMCDetector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {\r
- //\r
- // returns the momentum resolution \r
- //\r
- \r
- TGraph *graph = new TGraph(20, 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 * KMCDetector::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 ( 20, 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 ( 20, 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 * KMCDetector::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[20];\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
- KMCLayer *l = (KMCLayer*)fLayers.At(i);\r
- if (!l->IsDead() && l->fR>0) {\r
- layerRadius[count] = l->fR;\r
- layerThickness[count] = l->fx2X0;\r
- if (axis==0) {\r
- layerResolution[count] = l->fPhiRes;\r
- } else {\r
- layerResolution[count] = l->fZRes;\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 < 20 ; 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 ( 20, 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 * KMCDetector::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[20]; // with chosen particle mass\r
- Double_t kaonEfficiency[20], pionEfficiency[20], d0efficiency[20]; \r
- Double_t partEfficiency[2][20];\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 < 20 ; 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 < 20 ; j++ ) {\r
- d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);\r
- }\r
- } else { \r
- for ( Int_t j = 0 ; j < 20 ; 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 < 20 ; 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 ( 20, fTransMomenta, particleEfficiency ) ; // choosen mass\r
- graph->SetLineWidth(1);\r
- } else if (particle==1) {\r
- graph = new TGraph ( 20, fTransMomenta, pionEfficiency ) ;\r
- graph->SetLineWidth(1);\r
- } else if (particle ==2) {\r
- graph = new TGraph ( 20, fTransMomenta, kaonEfficiency ) ;\r
- graph->SetLineWidth(1);\r
- } else if (particle ==3) {\r
- graph = new TGraph ( 20, 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
-TGraph * KMCDetector::GetGraphRecoFakes(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
- //\r
-\r
- Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); \r
- \r
- Double_t particleFake[20]; // with chosen particle mass\r
- Double_t kaonFake[20], pionFake[20];\r
- Double_t partFake[2][20];\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 < 20 ; 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
- kaonFake[j] = fFake[1][j] /( doNotDecayFactor) ;\r
- } else { // PION\r
- pionFake[j] = fFake[0][j] ; \r
- }\r
- partFake[0][j] = pionFake[j];\r
- partFake[1][j] = kaonFake[j];\r
- } \r
- }\r
- \r
- } else { \r
- for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
- particleFake[j] = fFake[2][j];\r
- // NOTE: Decay factor (see kaon) should be included to be realiable\r
- }\r
- }\r
-\r
- for ( Int_t j = 0 ; j < 20 ; j++ ) { \r
- pionFake[j] *= 100;\r
- kaonFake[j] *= 100;\r
- particleFake[j] *= 100;\r
- }\r
- \r
- TGraph * graph = 0;\r
- if (particle==0) {\r
- graph = new TGraph ( 20, fTransMomenta, particleFake ) ; // choosen mass\r
- graph->SetLineWidth(1);\r
- } else if (particle==1) {\r
- graph = new TGraph ( 20, fTransMomenta, pionFake ) ;\r
- graph->SetLineWidth(1);\r
- } else if (particle ==2) {\r
- graph = new TGraph ( 20, fTransMomenta, kaonFake ) ;\r
- graph->SetLineWidth(1);\r
- } \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("Fake (%)") ;\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* KMCDetector::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* KMCDetector::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
- case 15:\r
- return GetGraphRecoFakes(0, color, linewidth); // Fake tracked particle\r
- case 16:\r
- return GetGraphRecoFakes(1, color, linewidth); // Fake pion\r
- case 17:\r
- return GetGraphRecoFakes(2, color, linewidth); // Fake kaon\r
- default:\r
- printf(" Error: chosen graph number not valid\n");\r
- }\r
- return 0;\r
-\r
-}\r
-\r
-void KMCDetector::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon, int setVer) {\r
- \r
- // All New configuration with X0 = 0.3 and resolution = 4 microns\r
- \r
- AddLayer((char*)"bpipe_its",2.0,0.0022, 0.092); // beam pipe, 0.5 mm Be\r
- AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation\r
- if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {\r
- printf("vertex will work as constraint: %.4f %.4f\n",fgVtxConstraint[0],fgVtxConstraint[1]);\r
- SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);\r
- }\r
- //\r
- // new ideal Pixel properties?\r
- Double_t x0 = 0.0050;\r
- Double_t resRPhi = 0.0006;\r
- Double_t resZ = 0.0006;\r
- Double_t xrho = 0.0116; // assume 0.5mm of silicon for eloss\r
- //\r
- if (flagMon) {\r
- x0 *= 3./5.;\r
- xrho *=3./5.;\r
- resRPhi = 0.0004;\r
- resZ = 0.0004;\r
- }\r
-\r
- double sclD = 1;\r
- double sclZ = 1;\r
-\r
- // default all new \r
- if (setVer<=0) {\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ); \r
- }\r
- else if (setVer==1) {\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 2.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 3.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its5", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- //\r
- /*\r
- UInt_t patt[] = {\r
- KMCTrackSummary::Bits(1,1,1),\r
- KMCTrackSummary::Bits(0,0,0,1,1),\r
- KMCTrackSummary::Bits(0,0,0,0,0,1,1)\r
- };\r
- RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\r
- */\r
- }\r
- else if (setVer==2) {\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its4", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its5", 33.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ); \r
- //\r
- /*\r
- UInt_t patt[] = {\r
- KMCTrackSummary::Bits(1,1,1,1),\r
- KMCTrackSummary::Bits(0,0,0,0,1,1),\r
- KMCTrackSummary::Bits(0,0,0,0,0,0,1,1,1)\r
- };\r
- RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\r
- */\r
- }\r
- /*\r
- // last 2 layers strips\r
- double resRPStr=0.0020, resZStr = 0.0830;\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6", 39.6 , x0, xrho, resRPStr, resZStr); \r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPStr, resZStr); \r
- */\r
- //*/\r
- /*\r
- // resolution scaled with R as res(2.2) * R/2.2\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ); \r
- KMCLayer* lr0 = 0;\r
- for (int i=0;i<GetNActiveITSLayers();i++) {\r
- KMCLayer *lr = GetActiveLayer(i);\r
- if (lr->IsVertex()) continue;\r
- if (lr0==0) {\r
- lr0=lr; \r
- // continue;\r
- }\r
- double scl = 5*lr->GetRadius()/lr0->GetRadius();\r
- SetResolution((char*)lr->GetName(), resRPhi*scl, resZ*scl*4);\r
- }\r
- */\r
- /*\r
- // 1st 2 layers double\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its3a", 6.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ); \r
- */\r
- /*\r
- // last 4 layers doubled\r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4a", 12.8 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5a", 23.9 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6a", 40.0 , x0, xrho, resRPhi, resZ); \r
-\r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its7a", 43.4 , x0, xrho, resRPhi, resZ); \r
- */\r
-\r
- /* //last 3 lr together \r
- AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its5", 42.2 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its6", 42.6 , x0, xrho, resRPhi, resZ); \r
- AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ); \r
- */\r
- if (flagTPC) {\r
- AddTPC(0.1,0.1); // TPC\r
- }\r
- PrintLayout();\r
-}\r
-\r
-void KMCDetector::MakeAliceCurrent(Bool_t flagTPC, Int_t AlignResiduals) {\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
- // 1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};\r
-\r
- AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation\r
- if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {\r
- printf("vertex will work as constraint: %.4f %.4f",fgVtxConstraint[0],fgVtxConstraint[1]);\r
- SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);\r
- }\r
- //\r
- AddLayer((char*)"bpipe_its",2.94,0.0022, 1.48e-01); // beam pipe\r
- AddLayer((char*)"tshld1_its",11.5,0.0065, 1.34e-01); // Thermal shield // 1.3% /2\r
- AddLayer((char*)"tshld2_its",31.0,0.0108, 2.22e-01); // Thermal shield // 1.3% /2\r
- //\r
- if (flagTPC) {\r
- AddTPC(0.1,0.1); // TPC\r
- }\r
- // Adding the ITS - current configuration\r
- \r
- if (AlignResiduals==0) {\r
- AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, 0.0012, 0.0130);\r
- AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, 0.0012, 0.0130);\r
- AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, 0.0035, 0.0025);\r
- AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, 0.0035, 0.0025);\r
- AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 0.0020, 0.0830);\r
- AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 0.0020, 0.0830);\r
- /*\r
- UInt_t patt[] = {\r
- KMCTrackSummary::Bits(1,1),\r
- KMCTrackSummary::Bits(0,0,1,1),\r
- KMCTrackSummary::Bits(0,0,0,0,1,1)\r
- };\r
- RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));\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_its", 3.9, 0.0114, 2.48e-01, 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_its", 7.6, 0.0114, 2.57e-01, 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_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),\r
- TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));\r
- AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 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_its",43.0, 0.0086, 2.25e-01, 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
- \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_its", 3.9, 0.0114, 2.48e-01, 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_ots", 7.6, 0.0114, 2.57e-01, 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_ots",15.0, 0.0113, 3.34e-01, 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_its",23.9, 0.0126, 3.50e-01, 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_its",38.0, 0.0083, 2.38e-01, 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_its",43.0, 0.0086, 2.25e-01, 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_its", 3.9, 0.0114, 2.48e-01, 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_its", 7.6, 0.0114, 2.57e-01, 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_its",15.0, 0.0113, 3.34e-01, 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_its",23.9, 0.0126, 3.50e-01, 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_its",38.0, 0.0083, 2.38e-01, 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_its",43.0, 0.0086, 2.25e-01, 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
-void KMCDetector::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
-void KMCDetector::ApplyMS(KMCProbe* trc, double x2X0) const\r
-{\r
- // simulate random modification of track params due to the MS\r
- if (x2X0<=0) return;\r
- double alpha = trc->GetAlpha(); // store original alpha\r
- double mass = trc->GetMass();\r
- //\r
- double snp = trc->GetSnp();\r
- double dip = trc->GetTgl();\r
- Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));\r
- x2X0 *= angle;\r
- //\r
- static double covCorr[15],covDum[21]={0};\r
- static double mom[3],pos[3];\r
- double *cov = (double*) trc->GetCovariance();\r
- memcpy(covCorr,cov,15*sizeof(double));\r
- trc->GetXYZ(pos);\r
- trc->GetPxPyPz(mom);\r
- double pt2 = mom[0]*mom[0]+mom[1]*mom[1];\r
- double pt = TMath::Sqrt(pt2);\r
- double ptot2 = pt2 + mom[2]*mom[2];\r
- double ptot = TMath::Sqrt(ptot2);\r
- double beta = ptot/TMath::Sqrt(ptot2 + mass*mass);\r
- double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);\r
- //\r
- // a la geant\r
- double phiSC = gRandom->Rndm()*TMath::Pi();\r
- double thtSC = gRandom->Gaus(0,1.4142*sigth);\r
- // printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);\r
- double sn = TMath::Sin(thtSC);\r
- double dx = sn*TMath::Sin(phiSC);\r
- double dy = sn*TMath::Cos(phiSC); \r
- double dz = TMath::Cos(thtSC);\r
- double v[3];\r
- // printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);\r
- for (int i=3;i--;) mom[i] /= ptot;\r
- double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);\r
- if (!IsZero(pt)) {\r
- double pd1 = mom[0]/vmm;\r
- double pd2 = mom[1]/vmm;\r
- v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;\r
- v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;\r
- v[2] = -vmm*dx + mom[2]*dz;\r
- }\r
- else {\r
- v[0] = dx;\r
- v[1] = dy;\r
- v[2] = dz*TMath::Sign(1.,mom[2]);\r
- }\r
- double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);\r
- // printf("before :%+e %+e %+e || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2], sigth, x2X0, pt, beta);\r
- // trc->Print();\r
- // direction cosines -> p\r
- for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;\r
- // printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);\r
- trc->Set(pos,mom,covDum,trc->Charge());\r
- //\r
- trc->Rotate(alpha);\r
- memcpy(cov,covCorr,15*sizeof(double));\r
- //\r
-}\r
-\r
-//________________________________________________________________________________\r
-Bool_t KMCDetector::TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const\r
-{\r
- // Transport track till layer maxLr, applying random MS\r
- //\r
- for (Int_t j=0; j<maxLr; j++) {\r
- KMCLayer* lr0 = (KMCLayer*)fLayers.At(j);\r
- KMCLayer* lr = (KMCLayer*)fLayers.At(j+1);\r
- if (!lr0->IsVertex()) {\r
- ApplyMS(probTr,lr0->fx2X0); // apply MS\r
- if (!probTr->CorrectForMeanMaterial(lr0,kFALSE)) return kFALSE; \r
- }\r
- //\r
- if (!PropagateToLayer(probTr,lr,1)) return kFALSE;\r
- // store randomized cluster local coordinates and phi\r
- lr->ResetBgClusters();\r
- double rz,ry;\r
- gRandom->Rannor(rz,ry);\r
- lr->GetMCCluster()->Set(probTr->GetY()+ry*lr->GetPhiRes(),probTr->GetZ()+rz*lr->GetZRes(), \r
- probTr->GetX(), probTr->GetAlpha() ); \r
- //\r
- }\r
- //\r
- return kTRUE;\r
-}\r
-\r
-//________________________________________________________________________________\r
-Bool_t KMCDetector::SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr,int nMC, int offset)\r
-{\r
- // analityc and fullMC (nMC trials) evaluaion of tracks with given kinematics.\r
- // the results are filled in KMCTrackSummary objects provided via summArr array\r
- //\r
- int progressP = 10;//0; // print progress in percents\r
- //\r
- progressP = int(nMC*0.01*progressP);\r
-\r
- if (!SolveSingleTrackViaKalman(mass,pt,eta)) return kFALSE;\r
- //\r
- // Store non-updated track errors of inward propagated seed >>>>>>>>\r
- int maxLr = fLastActiveITSLayer + offset;\r
- if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;\r
- KMCProbe probeTmp = fProbe; // original probe at vertex\r
- KMCLayer* lr = 0;\r
- for (Int_t j=1; j<=maxLr; j++) {\r
- lr = GetLayer(j);\r
- // printf("Here0: %d\n",j);\r
- if (!PropagateToLayer(&probeTmp,lr,1)) return 0;\r
- if (j!=maxLr) if (!probeTmp.CorrectForMeanMaterial(lr, kFALSE)) return 0;\r
- // printf("Prelim. Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(probeTmp.GetSigmaY2()),TMath::Sqrt(probeTmp.GetSigmaZ2()));\r
- }\r
- for (Int_t j=maxLr; j>0; j--) {\r
- lr = GetLayer(j);\r
- // printf("Here1: %d\n",j);\r
- if (j!=maxLr) if (!PropagateToLayer(&probeTmp,lr,-1)) return 0;\r
- lr->fSig2EstD = probeTmp.GetSigmaY2();\r
- lr->fSig2EstZ = probeTmp.GetSigmaZ2();\r
- // probeTmp.Print("l");\r
- printf("Natural Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(lr->fSig2EstD),TMath::Sqrt(lr->fSig2EstZ));\r
- if (!probeTmp.CorrectForMeanMaterial(lr, kTRUE)) return 0;\r
- }\r
- // Store non-updated track errors of inward propagated seed <<<<<<<<\r
- //\r
- int nsm = sumArr ? sumArr->GetEntriesFast() : 0;\r
- KMCLayer* vtx = GetLayer(0);\r
- //\r
- for (int i=0;i<nsm;i++) {\r
- KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(i);\r
- if (!tsm) continue;\r
- tsm->SetRefProbe( GetProbeTrack() ); // attach reference track (generated)\r
- tsm->SetAnProbe( vtx->GetAnProbe() ); // attach analitycal solution\r
- }\r
- //\r
- TStopwatch sw;\r
- sw.Start();\r
- for (int it=0;it<nMC;it++) {\r
- printf("ev: %d\n",it);\r
- SolveSingleTrackViaKalmanMC(offset);\r
- KMCProbe* trc = vtx->GetWinnerMCTrack();\r
- vtx->GetMCTracks()->Print();\r
- if (progressP==1 || (progressP>0 && (it%progressP)==0)) {\r
- printf("%d%% done |",it*100/nMC); \r
- sw.Stop(); sw.Print(); sw.Start(kFALSE);\r
- }\r
- for (int ism=nsm;ism--;) { // account the track in each of summaries\r
- KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(ism);\r
- if (!tsm) continue;\r
- tsm->AddUpdCalls(GetUpdCalls());\r
- tsm->AddTrack(trc); \r
- }\r
- }\r
- //\r
- sw.Stop();\r
- printf("Total time: "); sw.Print();\r
- return kTRUE;\r
-}\r
-\r
-//________________________________________________________________________________\r
-Int_t KMCDetector::GetLayerID(int actID) const\r
-{\r
- // find physical layer id from active id\r
- if (actID<0 || actID>fNActiveLayers) return -1;\r
- int start = actID<fNActiveITSLayers ? fLastActiveITSLayer : fLastActiveLayer;\r
- for (int i=start+1;i--;) {\r
- if (GetLayer(i)->GetActiveID()==actID) return i; \r
- }\r
- return -1;\r
-}\r
-\r
-//________________________________________________________________________________\r
-KMCProbe* KMCDetector::KalmanSmooth(int actLr, int actMin,int actMax) const\r
-{\r
- // estimate kalman smoothed track params at given active lr \r
- // from fit at layers actMin:actMax (excluding actLr)\r
- // SolveSingleTrackViaKalman must have been called before\r
- //\r
- if (actMin>actMax) swap(actMin,actMax);\r
- if (actMax>=fNActiveLayers) actMax = fNActiveLayers-1;\r
- int nlrfit = actMax-actMin;\r
- if (actLr>=actMin && actLr<=actMax) nlrfit-=1;\r
- if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}\r
- static KMCProbe iwd,owd;\r
- //\r
- // find phisical layer id's\r
- int pLr = GetLayerID(actLr);\r
- int pMin = GetLayerID(actMin);\r
- int pMax = GetLayerID(actMax);\r
- //\r
- // printf(">>> %d %d %d\n",pLr, pMin,pMax);\r
- Bool_t useIwd=kFALSE, useOwd=kFALSE;\r
- if (pLr<pMax) { // need inward piece\r
- iwd = GetLayer(pMax)->fTrCorr;\r
- iwd.ResetCovMat();\r
- iwd.GetHitsPatt() = 0;\r
- for (int i=pMax;i>=pLr;i--) {\r
- KMCLayer* lr = GetLayer(i);\r
- // printf("IWD %d\n",i);\r
- if (!lr->IsDead() && i!=pLr && i>=pMin) if (!UpdateTrack(&iwd,lr,&lr->fClCorr)) return 0;\r
- if (i!=pLr) {\r
- if (!iwd.CorrectForMeanMaterial(lr,kTRUE)) return 0; // correct for materials of this layer\r
- if (!PropagateToLayer(&iwd,GetLayer(i-1),-1)) return 0; // propagate to next layer\r
- }\r
- // printf("IWD%d: ",i); iwd.Print("l");\r
- }\r
- useIwd = kTRUE;\r
- }\r
- if (pLr>pMin) { // need outward piece\r
- owd = GetLayer(pMin)->fTrCorr;\r
- owd.ResetCovMat();\r
- owd.GetHitsPatt() = 0;\r
- for (int i=pMin;i<=pLr;i++) {\r
- KMCLayer* lr = GetLayer(i);\r
- // printf("OWD %d\n",i);\r
- if (!lr->IsDead() && i!=pLr && i<=pMax) if (!UpdateTrack(&owd,lr,&lr->fClCorr)) return 0;\r
- if (i!=pLr) {\r
- if (!owd.CorrectForMeanMaterial(lr,0)) return 0; // correct for materials of this layer\r
- if (!PropagateToLayer(&owd,GetLayer(i+1), 1)) return 0; // propagate to next layer\r
- }\r
- // printf("OWD%d: ",i); owd.Print("l");\r
- }\r
- useOwd = kTRUE;\r
- }\r
- //\r
- // was this extrapolation outside the fit range?\r
- if (!useIwd) return (KMCProbe*)&owd; \r
- if (!useOwd) return (KMCProbe*)&iwd;\r
- //\r
- // weight both tracks\r
- if (!iwd.Propagate(owd.GetAlpha(),owd.GetX(),fBFieldG)) return 0;\r
- double meas[2] = {owd.GetY(),owd.GetZ()};\r
- double measErr2[3] = {owd.GetSigmaY2(), owd.GetSigmaZY(), owd.GetSigmaZ2()};\r
- // printf("Weighting\n");\r
- // owd.Print("l");\r
- // iwd.Print("l");\r
- if (!iwd.Update(meas,measErr2)) return 0;\r
- iwd.GetHitsPatt() |= owd.GetHitsPatt();\r
-\r
- // printf("->\n");\r
- // iwd.Print("l");\r
-\r
- return (KMCProbe*)&iwd;\r
- //\r
-}\r
-\r
-//________________________________________________________________________________\r
-KMCProbe* KMCDetector::KalmanSmoothFull(int actLr, int actMin,int actMax) const\r
-{\r
- // estimate kalman smoothed track params at given active lr \r
- // from fit at layers actMin:actMax (excluding actLr)\r
- // SolveSingleTrackViaKalman must have been called before\r
- //\r
- static TClonesArray prediction("KMCProbe",10);\r
- static TClonesArray update("KMCProbe",10);\r
- static KMCProbe res;\r
- //\r
- if (actMin>actMax) swap(actMin,actMax);\r
- int nlrfit = actMax-actMin;\r
- if (actLr>=actMin && actLr<=actMax) nlrfit-=1;\r
- if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}\r
- //\r
- // find phisical layer id's\r
- int pLr = GetLayerID(actLr);\r
- int pMin = GetLayerID(actMin);\r
- int pMax = GetLayerID(actMax);\r
- //\r
- int dir=0,dirInt=0;\r
- if (pLr<=pMin) dir=-1; // inward extrapolation\r
- else if (pLr>=pMax) dir= 1; // outward extrapolation\r
- else if (actMax-actLr >= actLr-actMin) dirInt = -1; // inward interpolation (the test point is closer to inner layer)\r
- else dirInt = 1; // outward interpolation (the test point is closer to outer layer)\r
- //\r
- if (dir!=0) { // no sens to do smoothing: simple Kalman filtering extrapolation\r
- int start = dir<0 ? pMax : pMin;\r
- res = GetLayer(start)->fTrCorr;\r
- res.ResetCovMat();\r
- KMCLayer* lr = 0;\r
- for (int i=(dir<0?pMax:pMin); i!=pLr; i+=dir) { // track till nearest layer to pLr\r
- lr = GetLayer(i);\r
- if (!lr->IsDead() && !(i<pMin ||i>pMax)) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range\r
- if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer\r
- if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer \r
- }\r
- if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this nearest layer\r
- if (!PropagateToLayer(&res,GetLayer(pLr), dir)) return 0; // propagate to test layer\r
- return (KMCProbe*)&res;\r
- }\r
- //\r
- // too bad, need to do real filtering\r
- //\r
- int start = dirInt<0 ? pMax : pMin;\r
- int stop = dirInt<0 ? pMin-1 : pMax+1;\r
- res = GetLayer(start)->fTrCorr;\r
- res.ResetCovMat();\r
- KMCLayer* lr = 0;\r
- int count = 0;\r
- for (int i=start; i!=stop; i+=dirInt) { // track in full range, storing updates and predictions\r
- new(prediction[count]) KMCProbe(res);\r
- lr = GetLayer(i);\r
- if (!lr->IsDead() && i!=pLr) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range\r
- new(update[count]) KMCProbe(res);\r
- if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer\r
- if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer \r
- count++;\r
- }\r
- return (KMCProbe*)&res;\r
- //\r
-}\r
-\r
-//________________________________________________________________________________\r
-Bool_t KMCDetector::SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta)\r
-{\r
- // analytical estimate of tracking resolutions\r
- // fProbe.SetUseLogTermMS(kTRUE);\r
- //\r
- if (fMinITSHits>fNActiveITSLayers) {fMinITSHits = fNActiveITSLayers; printf("Redefined request of min N ITS hits to %d\n",fMinITSHits);}\r
- if (TMath::Abs(eta)<1e-3) fDensFactorEta = 1.;\r
- else {\r
- fDensFactorEta = TMath::Tan( 2.*TMath::ATan(TMath::Exp(-TMath::Abs(eta))) );\r
- fDensFactorEta = 1./TMath::Sqrt( 1. + 1./fDensFactorEta/fDensFactorEta);\r
- }\r
- double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta)); \r
- KMCProbe* probe = PrepareKalmanTrack(pt,lambda,mass,-1);\r
- if (!probe) return kFALSE;\r
- //\r
- KMCLayer *lr = 0;\r
- //\r
- //\r
- // Start the track fitting --------------------------------------------------------\r
- //\r
- // Back-propagate the covariance matrix along the track. \r
- // Kalman loop over the layers\r
- //\r
- KMCProbe* currTr = 0;\r
- lr = (KMCLayer*)fLayers.At(fLastActiveLayerTracked);\r
- lr->fTrCorr = *probe;\r
- delete probe; // rethink...\r
- //\r
- for (Int_t j=fLastActiveLayerTracked; j--; ) { // Layer loop\r
- //\r
- KMCLayer *lrP = lr;\r
- lr = (KMCLayer*)fLayers.At(j);\r
- //\r
- lr->fTrCorr = lrP->fTrCorr;\r
- currTr = &lr->fTrCorr;\r
- currTr->ResetHit(lrP->GetActiveID());\r
- //\r
- // if there was a measurement on prev layer, update the track\r
- if (!lrP->IsDead()) { // include measurement\r
- KMCCluster cl(currTr->GetY(),currTr->GetZ(), currTr->GetX(), currTr->GetAlpha());\r
- if (!UpdateTrack(currTr,lrP,&cl)) return kFALSE;\r
- }\r
- if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) return kFALSE; // correct for materials of this layer\r
- if (!PropagateToLayer(currTr,lr,-1)) return kFALSE; // propagate to current layer\r
- //\r
- } // end loop over layers\r
- //\r
- return kTRUE;\r
-}\r
-\r
-//____________________________________________________________\r
-Bool_t KMCDetector::SolveSingleTrackViaKalmanMC(int offset)\r
-{\r
- // MC estimate of tracking resolutions/effiencies. Requires that the SolveSingleTrackViaKalman\r
- // was called before, since it uses data filled by this method\r
- //\r
- // The MC tracking will be done starting from fLastActiveITSLayer + offset (before analytical estimate will be used)\r
- //\r
- // At this point, the fProbe contains the track params generated at vertex.\r
- // Clone it and propagate to target layer to generate hit positions affected by MS\r
- //\r
- fUpdCalls = 0.;\r
- KMCProbe *currTrP=0,*currTr=0;\r
- int maxLr = fLastActiveITSLayer + offset;\r
- if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;\r
- ResetMCTracks(maxLr);\r
- KMCLayer* lr = (KMCLayer*)fLayers.At(maxLr);\r
- currTr = lr->AddMCTrack(&fProbe); // start with original track at vertex\r
- //\r
- if (!TransportKalmanTrackWithMS(currTr, maxLr)) return kFALSE; // transport it to outermost layer where full MC is done\r
- //\r
- if (fLastActiveITSLayer<fLastActiveLayerTracked) { // prolongation from TPC\r
- // start from correct track propagated from above till maxLr\r
- double *covMS = (double*)currTr->GetCovariance();\r
- const double *covIdeal =lr->fTrCorr.GetCovariance();\r
- for (int i=15;i--;) covMS[i] = covIdeal[i];\r
- }\r
- else { // ITS SA: randomize the starting point\r
- // double *pars = (double*)currTr->GetParameter();\r
- // pars[0] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaY2()));\r
- // pars[1] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaZ2()));\r
- //\r
- currTr->ResetCovMat();\r
- /*\r
- double *trCov = (double*)currTr->GetCovariance();\r
- double *trPars = (double*)currTr->GetParameter();\r
- const double kLargeErr2PtI = 0.3*0.3;\r
- trCov[14] = TMath::Max(trCov[14],kLargeErr2PtI*trPars[4]*trPars[4]);\r
- */\r
- }\r
- //\r
- for (Int_t j=maxLr; j--; ) { // Layer loop\r
- //\r
- KMCLayer *lrP = lr;\r
- lr = (KMCLayer*)fLayers.At(j);\r
- int ntPrev = lrP->GetNMCTracks();\r
- //\r
- if (lrP->IsDead()) { // for passive layer just propagate the copy of all tracks of prev layer >>>\r
- for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
- currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
- currTr = lr->AddMCTrack( currTrP );\r
- if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer\r
- if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer\r
- }\r
- continue;\r
- } // treatment of dead layer <<<\r
- //\r
- if (lrP->IsTPC()) { // we don't consider bg hits in TPC, just update with MC cluster\r
- for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
- currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
- currTr = lr->AddMCTrack( currTrP );\r
- if (!UpdateTrack(currTr, lrP, lrP->GetMCCluster(), kTRUE)) {currTr->Kill(); continue;} // update with correct MC cl.\r
- if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer\r
- if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer\r
- }\r
- continue;\r
- } // treatment of ideal (TPC?) layer <<<\r
- //\r
- // active layer under eff. study (ITS?): propagate copy of every track to MC cluster frame (to have them all in the same frame)\r
- // and calculate the limits of bg generation\r
- KMCCluster* clMC = lrP->GetMCCluster();\r
- if (lrP->GetLayerEff()<gRandom->Rndm()) clMC->Kill(); // simulate inefficiency\r
- ResetSearchLimits();\r
- int nseeds = 0;\r
- for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer\r
- currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;\r
- currTr = lr->AddMCTrack( currTrP );\r
- if (!currTr->PropagateToCluster(clMC,fBFieldG)) {currTr->Kill(); continue;} // go to MC cluster\r
- if ( !(currTr->GetNITSHits()>0 && currTr->GetNITSHits()==currTr->GetNFakeITSHits()) ) UpdateSearchLimits(currTr, lrP); // RS\r
- nseeds++;\r
- }\r
- //\r
- // printf("%3d seeds\n",nseeds);\r
- if (fUseBackground && lrP->IsITS()) GenBgClusters(lrP); // generate background hits\r
- //\r
- ntPrev = lr->GetNMCTracks();\r
- for (int itr=ntPrev;itr--;) { // loop over all tracks PROPAGATED from previous layer to clusters frame on previous layer\r
- currTrP = lr->GetMCTrack(itr); // this is a seed from prev layer. The new clusters are attached to its copies, the seed itself\r
- // will be propagated w/o cluster update if it does not violate requested "reconstructed" track settings\r
- if (currTrP->IsKilled()) continue;\r
- //printf("Check %d %p %d\n",itr,currTrP,currTrP->GetUniqueID()); currTrP->Print();\r
- CheckTrackProlongations(currTrP, lr, lrP);\r
- if (NeedToKill(currTrP)) currTrP->Kill(); // kill track which was not updated at lrP\r
- //currTrP->Kill(); // kill track which was not updated at lrP\r
- }\r
- // \r
- lr->GetMCTracks()->Sort();\r
- int ntTot = lr->GetNMCTracks(); // propagate max amount of allowed tracks to current layer\r
- if (ntTot>fMaxSeedToPropagate && fMaxSeedToPropagate>0) {\r
- for (int itr=ntTot;itr>=fMaxSeedToPropagate;itr--) lr->GetMCTracks()->RemoveAt(itr);\r
- ntTot = fMaxSeedToPropagate;\r
- }\r
- //\r
- for (int itr=ntTot;itr--;) {\r
- currTr = lr->GetMCTrack(itr);\r
- if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill();continue;} // correct for materials of prev. layer\r
- if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill();continue;} // propagate to current layer\r
- }\r
- AliDebug(1,Form("Got %d tracks on layer %s",ntTot,lr->GetName()));\r
- // lr->GetMCTracks()->Print();\r
- //\r
- } // end loop over layers \r
- //\r
- // do we use vertex constraint?\r
- KMCLayer *vtx = GetLayer(0);\r
- if (!vtx->IsDead() && vtx->IsITS()) {\r
- int ntr = vtx->GetNMCTracks();\r
- for (int itr=0;itr<ntr;itr++) {\r
- currTr = vtx->GetMCTrack(itr);\r
- if (currTr->IsKilled()) continue;\r
- KMCCluster* clv = vtx->GetMCCluster();\r
- double meas[2] = {clv->GetY(),clv->GetZ()};\r
- double measErr2[3] = {vtx->fPhiRes*vtx->fPhiRes,0,vtx->fZRes*vtx->fZRes};\r
- double chi2v = currTr->GetPredictedChi2(meas,measErr2);\r
- currTr->AddHit(vtx->GetActiveID(), chi2v, -1);\r
- currTr->SetInnerLrChecked(vtx->GetActiveID());\r
- if (NeedToKill(currTr)) currTr->Kill();\r
- // if (vtx->IsITS()) {if (!UpdateTrack(currTr, vtx, vtx->GetMCCluster(), kFALSE)) {currTr->Kill();continue;}}\r
- }\r
- }\r
- EliminateUnrelated();\r
- \r
- return kTRUE;\r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t KMCDetector::PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const\r
-{\r
- // bring the track to layer and rotat to frame normal to its surface\r
- if (!trc->PropagateToR(lr->fR,fBFieldG, dir)) return kFALSE;\r
- //\r
- // rotate to frame with X axis normal to the surface (defined by ideal track)\r
- if (!lr->IsVertex()) {\r
- double pos[3];\r
- trc->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;\r
- if (!trc->Rotate(phi)) {\r
- AliDebug(2,Form("Failed to rotate to the frame (phi:%+.3f) of layer at %.2f at XYZ: %+.3f %+.3f %+.3f",\r
- phi,lr->fR,pos[0],pos[1],pos[2]));\r
- if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");\r
- return kFALSE;\r
- }\r
- }\r
- //\r
- return kTRUE;\r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t KMCDetector::UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster) const\r
-{\r
- // update track with measured cluster\r
- // propagate to cluster\r
- double meas[2] = {cl->GetY(),cl->GetZ()}; // ideal cluster coordinate\r
- double measErr2[3] = {lr->fPhiRes*lr->fPhiRes,0,lr->fZRes*lr->fZRes};\r
- //\r
- if (goToCluster) if (!trc->PropagateToCluster(cl,fBFieldG)) return kFALSE; // track was not propagated to cluster frame\r
- //\r
- double chi2 = trc->GetPredictedChi2(meas,measErr2);\r
- // if (chi2>fMaxChi2Cl) return kTRUE; // chi2 is too large\r
- // \r
- if (!trc->Update(meas,measErr2)) {\r
- AliDebug(2,Form("layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",\r
- lr->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));\r
- if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");\r
- return kFALSE;\r
- }\r
- trc->AddHit(lr->GetActiveID(), chi2);\r
- //\r
- return kTRUE;\r
-}\r
-\r
-//____________________________________________________________________________\r
-Int_t KMCDetector::GenBgClusters(KMCLayer* lr)\r
-{\r
- // Generate fake clusters in precalculated RPhi,Z range\r
- if (fNBgLimits<1) return 0; // limits were not set - no track was prolongated\r
- //\r
- // Fix search limits to avoid seeds which will anyway point very far from the vertex\r
- double tolY = TMath::Sqrt(lr->fSig2EstD)*fMaxChi2ClSQ;\r
- double tolZ = TMath::Sqrt(lr->fSig2EstZ)*fMaxChi2ClSQ;\r
- \r
- // printf("Before: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);\r
- if (fBgYMin < lr->fClCorr.fY-tolY) fBgYMin = lr->fClCorr.fY-tolY;\r
- if (fBgYMax > lr->fClCorr.fY+tolY) fBgYMax = lr->fClCorr.fY+tolY;\r
- if (fBgZMin < lr->fClCorr.fZ-tolZ) fBgZMin = lr->fClCorr.fZ-tolZ;\r
- if (fBgZMax > lr->fClCorr.fZ+tolZ) fBgZMax = lr->fClCorr.fZ+tolZ;\r
- //printf("After: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);\r
- //\r
- double dy = fBgYMax - fBgYMin;\r
- double dz = fBgZMax - fBgZMin;\r
- double surf = dy*dz; // surface of generation\r
- if (surf<0) return 0;\r
- double poissProb = surf*HitDensity(lr->fR)*lr->GetLayerEff();\r
- AliDebug(2,Form("Bg for Lr %s (r=%.2f) : Density %.2f on surface %.2e [%+.4f : %+.4f][%+.4f %+.4f]",\r
- lr->GetName(),lr->fR,HitDensity(lr->fR),surf,fBgYMin,fBgYMax,fBgZMin,fBgZMax));\r
- int nFakesGen = gRandom->Poisson( poissProb ); // preliminary number of extra clusters to test\r
- KMCCluster *refCl = lr->GetMCCluster();\r
- double sig2y = lr->GetPhiRes()*lr->GetPhiRes();\r
- double sig2z = lr->GetZRes()*lr->GetZRes();\r
- for (int ic=nFakesGen;ic--;) {\r
- double y = fBgYMin+dy*gRandom->Rndm();\r
- double z = fBgZMin+dz*gRandom->Rndm();\r
- double dfy = y-refCl->GetY();\r
- double dfz = z-refCl->GetZ();\r
- double dist = (dfy*dfy)/sig2y + (dfz*dfz)/sig2z;\r
- if (dist<4) continue; // avoid overlap with MC cluster\r
- lr->AddBgCluster(y, z, refCl->GetX(), refCl->GetPhi());\r
- }\r
- AliDebug(2,Form("Added %6d noise clusters on lr %s (poisson Prob=%8.2f for surface %.2e) DY:%7.4f DZ: %7.4f",\r
- lr->GetNBgClusters(),lr->GetName(),poissProb,surf,dy,dz));\r
- return nFakesGen;\r
- //\r
-}\r
-\r
-//____________________________________________________________________________\r
-void KMCDetector::UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr)\r
-{\r
- // define the search window for track on layer (where the bg hist will be generated)\r
- static double *currYMin = fBgYMinTr.GetArray();\r
- static double *currYMax = fBgYMaxTr.GetArray();\r
- static double *currZMin = fBgZMinTr.GetArray();\r
- static double *currZMax = fBgZMaxTr.GetArray();\r
- //\r
- double sizeY = probe->GetSigmaY2(), sizeZ = probe->GetSigmaZ2();\r
- // /*\r
- if (probe->GetNITSHits()<3) sizeY = 10*lr->fSig2EstD;\r
- if (probe->GetNITSHits()<2) sizeZ = 10*lr->fSig2EstZ;\r
- sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY+lr->fPhiRes*lr->fPhiRes); // max deviation in rphi to accept\r
- sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ+lr->fZRes*lr->fZRes); // max deviation in dz to accept\r
- // */\r
- //\r
- /*\r
- if (probe->GetNITSHits()<3) sizeY = 1./(1./lr->fSig2EstD + 1./sizeY);\r
- if (probe->GetNITSHits()<2) sizeZ = 1./(1./lr->fSig2EstZ + 1./sizeZ);\r
- sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY); // max deviation in rphi to accept\r
- sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ); // max deviation in dz to accept\r
- */\r
- //\r
- // if (sizeY>2) sizeY=2;\r
- // if (sizeZ>2) sizeZ=2;\r
- // printf("Sizes at %s: %.5f %.5f\n",lr->GetName(), sizeY,sizeZ);\r
- //\r
- if (fNBgLimits>=fBgYMinTr.GetSize()) { // expand arrays, update pointers\r
- fBgYMinTr.Set(2*(fNBgLimits+1));\r
- fBgYMaxTr.Set(2*(fNBgLimits+1));\r
- fBgZMinTr.Set(2*(fNBgLimits+1));\r
- fBgZMaxTr.Set(2*(fNBgLimits+1));\r
- currYMin = fBgYMinTr.GetArray();\r
- currYMax = fBgYMaxTr.GetArray();\r
- currZMin = fBgZMinTr.GetArray();\r
- currZMax = fBgZMaxTr.GetArray();\r
- }\r
- if (fBgYMin > (currYMin[fNBgLimits]=probe->GetY()-sizeY) ) fBgYMin = currYMin[fNBgLimits];\r
- if (fBgYMax < (currYMax[fNBgLimits]=probe->GetY()+sizeY) ) fBgYMax = currYMax[fNBgLimits];\r
- if (fBgZMin > (currZMin[fNBgLimits]=probe->GetZ()-sizeZ) ) fBgZMin = currZMin[fNBgLimits];\r
- if (fBgZMax < (currZMax[fNBgLimits]=probe->GetZ()+sizeZ) ) fBgZMax = currZMax[fNBgLimits];\r
- if (AliLog::GetGlobalDebugLevel()>=2) {\r
- probe->Print("l");\r
- AliInfo(Form("Seed%3d Lr %s limits for y:%+8.4f z:%+8.4f [%+.4f : %+.4f][%+.4f %+.4f]",fNBgLimits,lr->GetName(),probe->GetY(),probe->GetZ(),currYMin[fNBgLimits],currYMax[fNBgLimits],currZMin[fNBgLimits],currZMax[fNBgLimits]));\r
- AliInfo(Form("Global Limits Lr %s [%+.4f : %+.4f][%+.4f %+.4f]",lr->GetName(),fBgYMin,fBgYMax,fBgZMin,fBgZMax));\r
- AliInfo(Form("MC Cluster: %+.4f : %+.4f",lr->fClMC.fY, lr->fClMC.fZ));\r
- }\r
- probe->SetUniqueID(fNBgLimits++);\r
- //\r
- if (lr->IsITS() && probe->GetNFakeITSHits()==0) {\r
- if (fHMCLrResidRPhi) fHMCLrResidRPhi->Fill(probe->GetY() - lr->GetMCCluster()->GetY(), lr->GetActiveID());\r
- if (fHMCLrResidZ) fHMCLrResidZ->Fill(probe->GetZ() - lr->GetMCCluster()->GetZ(),lr->GetActiveID());\r
- }\r
- //\r
-}\r
-\r
-//____________________________________________________________________________\r
-void KMCDetector::CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP)\r
-{\r
- // explore prolongation of probe from lrP to lr with all possible clusters of lrP\r
- // the probe is already brought to clusters frame\r
- int nCl = lrP->GetNBgClusters();\r
- double measErr2[3] = {lrP->fPhiRes*lrP->fPhiRes,0,lrP->fZRes*lrP->fZRes};\r
- double meas[2] = {0,0};\r
- UInt_t tmpID = probe->GetUniqueID();\r
- double yMin = fBgYMinTr[tmpID];\r
- double yMax = fBgYMaxTr[tmpID];\r
- double zMin = fBgZMinTr[tmpID];\r
- double zMax = fBgZMaxTr[tmpID];\r
- //\r
- probe->SetInnerLrChecked(lrP->GetActiveID());\r
- for (int icl=-1;icl<nCl;icl++) {\r
- KMCCluster* cl = icl<0 ? lrP->GetMCCluster() : lrP->GetBgCluster(icl); // -1 is for true MC cluster\r
- if (cl->IsKilled()) {\r
- if (AliLog::GetGlobalDebugLevel()>1) {printf("Skip cluster %d ",icl); cl->Print();}\r
- continue;\r
- }\r
- double y = cl->GetY();\r
- double z = cl->GetZ();\r
- AliDebug(2,Form("Check seed%d against cl#%d out of %d at layer %s | y:%+8.4f z:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f]",tmpID,icl,nCl,lrP->GetName(),y,z,yMin,yMax,zMin,zMax));\r
- if (AliLog::GetGlobalDebugLevel()>0) {\r
- if (icl==-1 && probe->GetNFakeITSHits()==0) {\r
- meas[0] = y; meas[1] = z;\r
- double chi2a = probe->GetPredictedChi2(meas,measErr2);\r
- if (chi2a>fMaxChi2Cl || (y<yMin || y>yMax) || (z<zMin || z>zMax)) {\r
- probe->Print();\r
- printf("Loosing good point (y:%+8.4f z:%+8.4f) on lr %s: chi2: %.2f | dy:%+8.4f dz:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f] |x: %.2f %.2f | phi: %.2f %.2f\n",\r
- y,z,lrP->GetName(),chi2a,y-probe->GetY(),z-probe->GetZ(),yMin,yMax,zMin,zMax, probe->GetX(), cl->GetX(), probe->GetAlpha(), cl->GetPhi());\r
- }\r
- }\r
- }\r
- if (y<yMin || y>yMax) continue; // preliminary check on Y\r
- if (z<zMin || z>zMax) continue; // preliminary check on Z\r
- meas[0] = y; meas[1] = z;\r
- double chi2 = probe->GetPredictedChi2(meas,measErr2);\r
- if (fHMCLrChi2 && probe->GetNFakeITSHits()==0 && icl==-1) fHMCLrChi2->Fill(chi2,lrP->GetActiveID());\r
- AliDebug(2,Form("Seed-to-cluster chi2 = Chi2=%.2f",chi2));\r
- if (chi2>fMaxChi2Cl) continue;\r
- // \r
- // update track copy\r
- KMCProbe* newTr = lr->AddMCTrack( probe );\r
- fUpdCalls++;\r
- if (!newTr->Update(meas,measErr2)) {\r
- AliDebug(2,Form("Layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",\r
- lrP->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));\r
- if (AliLog::GetGlobalDebugLevel()>1) newTr->Print("l");\r
- newTr->Kill();\r
- continue;\r
- }\r
- newTr->AddHit(lrP->GetActiveID(), chi2, icl);\r
- if (AliLog::GetGlobalDebugLevel()>1) {\r
- AliInfo("Cloned updated track is:");\r
- newTr->Print();\r
- }\r
- if (NeedToKill(newTr)) newTr->Kill();\r
- }\r
- //\r
-}\r
-\r
-//____________________________________________________________________________\r
-void KMCDetector::ResetMCTracks(Int_t maxLr)\r
-{\r
- int nl = GetNLayers();\r
- if (maxLr<0 || maxLr>=nl) maxLr = nl-1;\r
- for (int i=maxLr+1;i--;) GetLayer(i)->ResetMCTracks();\r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t KMCDetector::NeedToKill(KMCProbe* probe) const\r
-{\r
- // check if the seed at given layer (last one where update was tried) \r
- // still has chances to be reconstructed\r
- const Bool_t kModeKillMiss = kFALSE;\r
- //\r
- Bool_t kill = kFALSE;\r
- while (1) {\r
- int il = probe->GetInnerLayerChecked();\r
- int nITS = probe->GetNITSHits();\r
- int nITSMax = nITS + il; // maximum it can have\r
- if (nITSMax<fMinITSHits) {\r
- kill = kTRUE; \r
- break;\r
- } // has no chance to collect enough ITS hits\r
- //\r
- int ngr = fPattITS.GetSize();\r
- if (ngr>0) { // check pattern\r
- UInt_t patt = probe->GetHitsPatt();\r
- // complete the layers not checked yet\r
- for (int i=il;i--;) patt |= (0x1<<i);\r
- for (int ig=ngr;ig--;) \r
- if (!(((UInt_t)fPattITS[ig]) & patt)) {\r
- kill = kTRUE; \r
- break;\r
- }\r
- //\r
- }\r
- //\r
- if (nITS>2) { // check if smallest possible norm chi2/ndf is acceptable\r
- double chi2min = probe->GetChi2();\r
- if (kModeKillMiss) {\r
- int nMiss = fNActiveITSLayers - probe->GetInnerLayerChecked() - nITS; // layers already missed\r
- chi2min = nMiss*probe->GetMissingHitPenalty();\r
- }\r
- chi2min /= ((nITSMax<<1)-KMCProbe::kNDOF);\r
- if (chi2min>fMaxNormChi2NDF) {\r
- kill = kTRUE; \r
- break;\r
- }\r
- }\r
- //\r
- // loose vertex constraint\r
- double dst;\r
- if (nITS>=2) {\r
- probe->GetZAt(0,fBFieldG,dst);\r
- //printf("Zd (F%d): %f\n",probe->GetNFakeITSHits(),dst);\r
- if (TMath::Abs(dst)>10.) {\r
- kill = kTRUE; \r
- break;\r
- }\r
- }\r
- if (nITS>=3) {\r
- probe->GetYAt(0,fBFieldG,dst);\r
- //printf("Dd (F%d): %f\n",probe->GetNFakeITSHits(),dst);\r
- if (TMath::Abs(dst)>10.) {\r
- kill = kTRUE; \r
- break;\r
- }\r
- }\r
- //\r
- break;\r
- }\r
- if (kill && AliLog::GetGlobalDebugLevel()>1 && probe->GetNFakeITSHits()==0) {\r
- printf("Killing good seed, last upd layer was %d\n",probe->GetInnerLayerChecked());\r
- probe->Print("l");\r
- }\r
- return kill;\r
-}\r
-\r
-//_____________________________________________________________________\r
-void KMCDetector::EliminateUnrelated()\r
-{\r
- // kill useless tracks\r
- KMCLayer* lr = GetLayer(0);\r
- int ntr = lr->GetNMCTracks();\r
- int nval = 0;\r
- for (int itr=0;itr<ntr;itr++) {\r
- KMCProbe* probe = lr->GetMCTrack(itr);\r
- if (probe->IsKilled()) continue;\r
- if (probe->GetNITSHits()-probe->GetNFakeITSHits()<1) {probe->Kill(); continue;}\r
- nval++;\r
- }\r
- lr->GetMCTracks()->Sort();\r
- const int kDump = 0;\r
- if (kDump>0) {\r
- printf("Valid %d out of %d\n",nval, ntr);\r
- ntr = ntr>kDump ? kDump:0;\r
- for (int itr=0;itr<ntr;itr++) {\r
- lr->GetMCTrack(itr)->Print();\r
- }\r
- }\r
-}\r
-\r
-//_____________________________________________________________________\r
-void KMCDetector::RequirePattern(UInt_t *patt, int groups)\r
-{\r
- if (groups<1) {fPattITS.Set(0); return;}\r
- fPattITS.Set(groups);\r
- for (int i=0;i<groups;i++) fPattITS[i] = patt[i];\r
-}\r
-\r
-//_____________________________________________________________________\r
-void KMCDetector::CalcHardSearchLimits(double dzv)\r
-{\r
- //\r
- TArrayD zlims;\r
- zlims.Set(fNActiveITSLayers);\r
- for (int il0=0;il0<fNActiveITSLayers;il0++) {\r
- KMCLayer* lr0 = GetActiveLayer(il0);\r
- double angZTol = dzv/lr0->GetRadius();\r
- for (int il1=0;il1<fNActiveITSLayers;il1++) {\r
- if (il1==il0) continue;\r
- KMCLayer* lr1 = GetActiveLayer(il1);\r
- double ztol = angZTol*TMath::Abs(lr0->GetRadius() - lr1->GetRadius());\r
- if (ztol>zlims[il1]) zlims[il1] = ztol;\r
- }\r
- }\r
- //\r
- for (int il=0;il<fNActiveITSLayers;il++) printf("ZTol%d: %8.4f\n",il,zlims[il]);\r
-}\r
-\r
-//_______________________________________________________\r
-double KMCDetector::PropagateBack(KMCProbe* trc) \r
-{\r
- static KMCProbe bwd;\r
- bwd = *trc;\r
- bwd.ResetCovMat();\r
- static double measErr2[3] = {0,0,0};\r
- static double meas[2] = {0,0};\r
- int icl = 0;\r
- double chi2Tot = 0;\r
- for (int il=1;il<=fLastActiveITSLayer;il++) {\r
- KMCLayer* lr = GetLayer(il);\r
- if (!PropagateToLayer(&bwd,lr,1)) return -1;\r
- int aID = lr->GetActiveID();\r
- if (aID>-1 && (icl=bwd.fClID[aID])>=-1) {\r
- KMCCluster* clMC = icl<0 ? lr->GetMCCluster() : lr->GetBgCluster(icl);\r
- if (!bwd.PropagateToCluster(clMC,fBFieldG)) return -1;\r
- meas[0] = clMC->GetY(); meas[1] = clMC->GetZ();\r
- measErr2[0] = lr->fPhiRes*lr->fPhiRes;\r
- measErr2[2] = lr->fZRes*lr->fZRes;\r
- double chi2a = bwd.GetPredictedChi2(meas,measErr2);\r
- chi2Tot += chi2a;\r
- printf("Chis %d (cl%+3d): t2c: %6.3f tot: %6.3f\n",aID,icl,chi2a, chi2Tot);\r
- bwd.Update(meas,measErr2);\r
- bwd.AddHit(aID, chi2a, icl); \r
- }\r
- if (!bwd.CorrectForMeanMaterial(lr,kFALSE)) return -1;\r
- }\r
- return chi2Tot;\r
-}\r
-\r
-\r
-//////////////////////////////////////////////////////////////////////////////////////////////\r
-//--------------------------------------------------------------------------------------------\r
-ClassImp(KMCTrackSummary)\r
-\r
-Int_t KMCTrackSummary::fgSumCounter = 0;\r
-\r
-//____________________________________________________________________________\r
-KMCTrackSummary::KMCTrackSummary(const char* name,const char* title, int nlr) :\r
- TNamed(name,name), fNITSLayers(nlr),\r
- fPattITS(0),fPattITSFakeExcl(0),fPattITSCorr(0),\r
- fHMCChi2(0),fHMCSigDCARPhi(0),fHMCSigDCAZ(0),fHMCSigPt(0),fHMCNCl(0),fHMCFakePatt(0),\r
- fCountAcc(0),fCountTot(0),fUpdCalls(0),\r
- fRefProbe(),fAnProbe()\r
-{\r
- // create summary structure for specific track conditions\r
- //\r
- SetMinMaxClITS();\r
- SetMinMaxClITSFake();\r
- SetMinMaxClITSCorr();\r
- //\r
- if (name==0 && title==0 && nlr==0) return; // default dummy contructor\r
- //\r
- enum {kNBinRes=1000};\r
- const double kMaxChi2(10),kMaxResRPhi=0.1, kMaxResZ=0.1,kMaxResPt=0.3;\r
- //\r
- TString strN = name, strT = title;\r
- if (strN.IsNull()) strN = Form("TrSum%d",fgSumCounter);\r
- if (strT.IsNull()) strT = strN;\r
- fgSumCounter++;\r
- //\r
- if (fNITSLayers<1) {\r
- fNITSLayers=KMCProbe::GetNITSLayers(); \r
- if (fNITSLayers<1) {AliError("N ITS layer is not provided and not available from KMCProbe::GetNITSLayers()"); exit(1);}\r
- AliInfo(Form("%s No N Layers provided, set %d",strN.Data(),fNITSLayers));\r
- }\r
- nlr = fNITSLayers;\r
- //\r
- TString nam,tit;\r
- //\r
- nam = Form("%s_mc_chi2",strN.Data());\r
- tit = Form("%s MC #chi^{2}",strT.Data());\r
- fHMCChi2 = new TH1F(nam.Data(),tit.Data(),kNBinRes,0,kMaxChi2);\r
- fHMCChi2->GetXaxis()->SetTitle("#chi^{2}");\r
- fHMCChi2->Sumw2();\r
- //\r
- nam = Form("%s_mc_DCArphi",strN.Data());\r
- tit = Form("%s MC #DeltaDCA r#phi",strT.Data());\r
- fHMCSigDCARPhi = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResRPhi,kMaxResRPhi);\r
- fHMCSigDCARPhi->GetXaxis()->SetTitle("#Delta r#phi");\r
- fHMCSigDCARPhi->Sumw2();\r
- //\r
- nam = Form("%s_mc_DCAz",strN.Data());\r
- tit = Form("%s MC #DeltaDCA Z",strT.Data());\r
- fHMCSigDCAZ = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResZ,kMaxResZ);\r
- fHMCSigDCAZ->GetXaxis()->SetTitle("#Delta Z");\r
- fHMCSigDCAZ->Sumw2();\r
- //\r
- nam = Form("%s_mc_pt",strN.Data());\r
- tit = Form("%s MC $Deltap_{T}/p_{T}",strT.Data());\r
- fHMCSigPt = new TH1F(nam.Data(),tit.Data(),2*kNBinRes,-kMaxResPt,kMaxResPt);\r
- fHMCSigPt->GetXaxis()->SetTitle("#Deltap_{T}/p_{T}");\r
- fHMCSigPt->Sumw2();\r
- //\r
- nam = Form("%s_n_cl",strN.Data());\r
- tit = Form("%s N Clusters",strT.Data());\r
- fHMCNCl = new TH2F(nam.Data(),tit.Data(),nlr,1,nlr+1,nlr,1,nlr+1);\r
- fHMCNCl->GetXaxis()->SetTitle("N Clusters Tot");\r
- fHMCNCl->GetYaxis()->SetTitle("N Clusters Fake");\r
- fHMCNCl->Sumw2();\r
- //\r
- nam = Form("%s_fake_cl",strN.Data());\r
- tit = Form("%s Fake Clusters",strT.Data());\r
- fHMCFakePatt = new TH1F(nam.Data(),tit.Data(),nlr,-0.5,nlr-0.5);\r
- fHMCFakePatt->GetXaxis()->SetTitle("Fake clusters pattern");\r
- fHMCFakePatt->Sumw2();\r
- //\r
-}\r
-\r
-//____________________________________________________________________________\r
-KMCTrackSummary::~KMCTrackSummary()\r
-{\r
- if (fHMCChi2) delete fHMCChi2;\r
- if (fHMCSigDCARPhi) delete fHMCSigDCARPhi;\r
- if (fHMCSigDCAZ) delete fHMCSigDCAZ;\r
- if (fHMCSigPt) delete fHMCSigPt;\r
- if (fHMCNCl) delete fHMCNCl;\r
- if (fHMCFakePatt) delete fHMCFakePatt;\r
-}\r
-\r
-\r
-//____________________________________________________________________________\r
-void KMCTrackSummary::Add(const KMCTrackSummary* src, double scl)\r
-{\r
- if (fHMCChi2) fHMCChi2->Add(src->fHMCChi2,scl);\r
- if (fHMCSigDCARPhi) fHMCSigDCARPhi->Add(src->fHMCSigDCARPhi,scl);\r
- if (fHMCSigDCAZ) fHMCSigDCAZ->Add(src->fHMCSigDCAZ,scl);\r
- if (fHMCSigPt) fHMCSigPt->Add(src->fHMCSigPt,scl);\r
- if (fHMCNCl) fHMCNCl->Add(src->fHMCNCl,scl);\r
- if (fHMCFakePatt) fHMCFakePatt->Add(src->fHMCFakePatt,scl);\r
- fCountAcc += src->fCountAcc;\r
- fUpdCalls += src->fUpdCalls;\r
-}\r
-\r
-//____________________________________________________________________________\r
-void KMCTrackSummary::PrependTNamed(TNamed* obj, const char* nm, const char* tit)\r
-{\r
- // prepend name, title by this prefix\r
- TString name = obj->GetName(),title = obj->GetTitle();\r
- name.Prepend(nm); title.Prepend(tit);\r
- obj->SetNameTitle(name.Data(),title.Data());\r
- //\r
-}\r
- \r
-//____________________________________________________________________________\r
-void KMCTrackSummary::SetNamePrefix(const char* pref)\r
-{\r
- // prepend all names by this prefix\r
- TString nmT,nmP = pref;\r
- if (nmP.IsNull()) return;\r
- nmP = Form("%s_",pref);\r
- nmT = Form("%s ",pref);\r
- PrependTNamed(this, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCChi2, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCSigDCARPhi, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCSigDCAZ, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCSigPt, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCNCl, nmP.Data(), nmT.Data());\r
- PrependTNamed(fHMCFakePatt, nmP.Data(), nmT.Data());\r
- // \r
-}\r
-\r
-//____________________________________________________________________________\r
-Bool_t KMCTrackSummary::CheckTrack(KMCProbe* trc)\r
-{\r
- // check if the track satisfies to selection conditions\r
- fCountTot++;\r
- if (!trc) return kFALSE;\r
- if (OutOfRange(trc->GetNITSHits(), fMinMaxClITS)) return kFALSE;\r
- if (OutOfRange(trc->GetNFakeITSHits(), fMinMaxClITSFake)) return kFALSE;\r
- if (OutOfRange(trc->GetNITSHits()-trc->GetNFakeITSHits(), fMinMaxClITSCorr)) return kFALSE;\r
- //\r
- // check layer patterns if requested\r
- UInt_t patt = trc->GetHitsPatt();\r
- UInt_t pattF = trc->GetFakesPatt();\r
- UInt_t pattC = patt&(~pattF); // correct hits\r
- // is there at least one hit in each of requested layer groups?\r
- for (int ip=fPattITS.GetSize();ip--;) if (!CheckPattern(patt,fPattITS[ip])) return kFALSE;\r
- // is there at least one fake in any of requested layer groups?\r
- for (int ip=fPattITSFakeExcl.GetSize();ip--;) if (CheckPattern(pattF,fPattITSFakeExcl[ip])) return kFALSE;\r
- // is there at least one correct hit in each of requested layer groups?\r
- for (int ip=fPattITSCorr.GetSize();ip--;) if (!CheckPattern(pattC,fPattITSCorr[ip])) return kFALSE;\r
- //\r
- fCountAcc++;\r
- return kTRUE;\r
-}\r
-\r
-//____________________________________________________________________________\r
-void KMCTrackSummary::AddTrack(KMCProbe* trc)\r
-{\r
- // fill track info\r
- if (!CheckTrack(trc)) return;\r
- //\r
- fHMCChi2->Fill(trc->GetNormChi2(kFALSE));\r
- fHMCSigDCARPhi->Fill(trc->GetY());\r
- fHMCSigDCAZ->Fill(trc->GetZ());\r
- if (fRefProbe.Pt()>0) fHMCSigPt->Fill( trc->Pt()/fRefProbe.Pt()-1.);\r
- // printf("Pts: %.3f %.3f -> %.3f\n",trc->Pt(),fRefProbe.Pt(),trc->Pt()/fRefProbe.Pt()-1.);\r
- // trc->Print("l");\r
- // fRefProbe.Print("l");\r
- fHMCNCl->Fill(trc->GetNITSHits(),trc->GetNFakeITSHits());\r
- for (int i=trc->GetNITSLayers();i--;) if (trc->IsHitFake(i)) fHMCFakePatt->Fill(i);\r
- //\r
-}\r
-\r
-//____________________________________________________________________________\r
-UInt_t KMCTrackSummary::Bits(Bool_t l0, Bool_t l1, Bool_t l2, Bool_t l3, Bool_t l4, Bool_t l5, Bool_t l6, Bool_t l7,Bool_t l8, Bool_t l9,\r
- Bool_t l10,Bool_t l11,Bool_t l12,Bool_t l13,Bool_t l14,Bool_t l15,Bool_t l16,Bool_t l17,Bool_t l18,Bool_t l19,\r
- Bool_t l20,Bool_t l21,Bool_t l22,Bool_t l23,Bool_t l24,Bool_t l25,Bool_t l26,Bool_t l27,Bool_t l28,Bool_t l29,\r
- Bool_t l30,Bool_t l31)\r
-{\r
- // create corresponding bit pattern\r
- UInt_t patt = 0;\r
- if (l0 ) patt |= (0x1<<0); if (l1 ) patt |= (0x1<<1); if (l2 ) patt |= (0x1<<2);\r
- if (l3 ) patt |= (0x1<<3); if (l4 ) patt |= (0x1<<4); if (l5 ) patt |= (0x1<<5);\r
- if (l6 ) patt |= (0x1<<6); if (l7 ) patt |= (0x1<<7); if (l8 ) patt |= (0x1<<8);\r
- if (l9 ) patt |= (0x1<<9); if (l10) patt |= (0x1<<10); if (l11) patt |= (0x1<<11);\r
- if (l12) patt |= (0x1<<12); if (l13) patt |= (0x1<<13); if (l14) patt |= (0x1<<14);\r
- if (l15) patt |= (0x1<<15); if (l16) patt |= (0x1<<16); if (l17) patt |= (0x1<<17);\r
- if (l18) patt |= (0x1<<18); if (l19) patt |= (0x1<<19); if (l20) patt |= (0x1<<20);\r
- if (l21) patt |= (0x1<<21); if (l22) patt |= (0x1<<22); if (l23) patt |= (0x1<<23);\r
- if (l24) patt |= (0x1<<24); if (l25) patt |= (0x1<<25); if (l26) patt |= (0x1<<26);\r
- if (l27) patt |= (0x1<<27); if (l28) patt |= (0x1<<28); if (l29) patt |= (0x1<<29);\r
- if (l30) patt |= (0x1<<30); if (l31) patt |= (0x1<<31); \r
- return patt;\r
-}\r
-\r
-//__________________________________________________________________________\r
-void KMCTrackSummary::Print(Option_t* ) const\r
-{\r
- printf("%s: summary for track M=%5.3f pT: %6.3f eta: %.2f\n", GetName(),\r
- fRefProbe.GetMass(),fRefProbe.Pt(), fRefProbe.Eta());\r
- printf("Cuts on NCl ITS: Tot: %2d - %2d Fake: %2d - %2d Corr: %2d - %2d\n",\r
- fMinMaxClITS[0],fMinMaxClITS[1], \r
- fMinMaxClITSFake[0],fMinMaxClITSFake[1],\r
- fMinMaxClITSCorr[0],fMinMaxClITSCorr[1]);\r
- //\r
- int nlr = fNITSLayers;\r
- if (fPattITS.GetSize()) {\r
- printf("Require at least 1 hit in groups: ");\r
- printf("Hits obligatory in groups: ");\r
- for (int i=fPattITS.GetSize();i--;) {\r
- UInt_t pat = (UInt_t)fPattITS[i];\r
- printf("[");\r
- for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
- printf("] ");\r
- }\r
- printf("\n");\r
- }\r
- //\r
- if (fPattITSFakeExcl.GetSize()) {\r
- printf("Fake hit forbidden in groups : ");\r
- for (int i=fPattITSFakeExcl.GetSize();i--;) {\r
- UInt_t pat = (UInt_t)fPattITSFakeExcl[i];\r
- printf("[");\r
- for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
- printf("] ");\r
- }\r
- printf("\n");\r
- }\r
- //\r
- if (fPattITSCorr.GetSize()) {\r
- printf("Correct hit obligatory in groups: ");\r
- for (int i=fPattITSCorr.GetSize();i--;) {\r
- UInt_t pat = (UInt_t)fPattITSCorr[i];\r
- printf("[");\r
- for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);\r
- printf("] ");\r
- }\r
- printf("\n");\r
- }\r
- // \r
- printf("Entries: Tot: %.4e Acc: %.4e -> Eff: %.4f+-%.4f %.2e KMC updates\n",fCountTot,fCountAcc,GetEff(),GetEffErr(),GetUpdCalls());\r
-}\r
-\r
+#include "KMCDetector.h"
+#include <TMath.h>
+#include <TH1F.h>
+#include <TProfile.h>
+#include <TF1.h>
+#include <TMatrixD.h>
+#include <TGraph.h>
+#include <TAxis.h>
+#include <TFormula.h>
+#include <TCanvas.h>
+#include <TEllipse.h>
+#include <TText.h>
+#include <TRandom3.h>
+#include <TGraphErrors.h>
+#include <TStopwatch.h>
+
+/***********************************************************
+
+Fast Simulation tool for Inner Tracker Systems
+
+***********************************************************/
+
+
+#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
+
+ClassImp(KMCProbe)
+
+Int_t KMCProbe::fgNITSLayers = 0;
+Double_t KMCProbe::fgMissingHitPenalty = 2.;
+//__________________________________________________________________________
+KMCProbe& KMCProbe::operator=(const KMCProbe& src)
+{
+ if (this!=&src) {
+ AliExternalTrackParam::operator=(src);
+ fMass = src.fMass;
+ fChi2 = src.fChi2;
+ fHits = src.fHits;
+ fFakes = src.fFakes;
+ fNHits = src.fNHits;
+ fNHitsITS = src.fNHitsITS;
+ fNHitsITSFake = src.fNHitsITSFake;
+ fInnLrCheck = src.fInnLrCheck;
+ for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
+ }
+ return *this;
+}
+
+//__________________________________________________________________________
+KMCProbe::KMCProbe(KMCProbe& src)
+ : AliExternalTrackParam(src),
+ fMass(src.fMass),
+ fChi2(src.fChi2),
+ fHits(src.fHits),
+ fFakes(src.fFakes),
+ fNHits(src.fNHits),
+ fNHitsITS(src.fNHitsITS),
+ fNHitsITSFake(src.fNHitsITSFake),
+ fInnLrCheck(src.fInnLrCheck)
+{
+ for (int i=kMaxITSLr;i--;) fClID[i] = src.fClID[i];
+}
+
+//__________________________________________________________________________
+void KMCProbe::ResetCovMat()
+{
+ // reset errors
+ double *trCov = (double*)GetCovariance();
+ double *trPars = (double*)GetParameter();
+ const double kLargeErr2Coord = 50*50;
+ const double kLargeErr2Dir = 0.6*0.6;
+ const double kLargeErr2PtI = 0.5*0.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];
+ //
+}
+
+//__________________________________________________________________________
+void KMCProbe::Print(Option_t* option) const
+{
+ printf("M=%.3f Chi2=%7.2f (Norm:%6.2f|%d) Hits: Total:%d ITS:%d ITSFakes:%d | Y:%+8.4f Z: %+8.4f |",
+ fMass,fChi2,GetNormChi2(kTRUE),fInnLrCheck,fNHits,fNHitsITS,fNHitsITSFake, GetY(),GetZ());
+ for (int i=0;i<fgNITSLayers;i++) {
+ if (!IsHit(i)) printf(".");
+ else printf("%c",IsHitFake(i) ? '-':'+');
+ }
+ printf("|%s\n",IsKilled() ? " KILLED":"");
+ TString opt = option;
+ if (!opt.IsNull()) AliExternalTrackParam::Print(option);
+}
+
+//__________________________________________________________________________
+Int_t KMCProbe::Compare(const TObject* obj) const
+{
+ // compare to tracks
+ const KMCProbe* trc = (KMCProbe*) obj;
+ if (trc->IsKilled()) {
+ if (IsKilled()) return 0;
+ return -1;
+ }
+ else if (IsKilled()) return 1;
+ double chi2a = GetNormChi2(kTRUE);
+ double chi2b = trc->GetNormChi2(kTRUE);
+ if (chi2a<chi2b) return -1;
+ if (chi2a>chi2b) return 1;
+ return 0;
+}
+
+//__________________________________________________________________________
+Bool_t KMCProbe::GetXatLabR(Double_t r,Double_t &x, Double_t bz, Int_t dir) const
+{
+ // 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 = GetParameter();
+ const Double_t &fy=pars[0], &sn = pars[2];
+ //
+ double fx = GetX();
+ double crv = 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\n",x0,y0);
+ //
+ if (r0<=kAlmost0) {
+ AliDebug(2,Form("r0 = %f",r0));
+ 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) {
+ AliDebug(2,Form("g=%f tR=%f r0=%f\n",g,tR, r0));
+ 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;
+}
+
+//____________________________________
+Bool_t KMCProbe::PropagateToR(double r, double b, int dir)
+{
+ // go to radius R
+ //
+ double xR = 0;
+ double rr = r*r;
+ int iter = 0;
+ const double kTiny = 1e-4;
+ while(1) {
+ if (!GetXatLabR(r ,xR, b, dir)) {
+ // printf("Track with pt=%f cannot reach radius %f\n",Pt(),r);
+ // Print("l");
+ return kFALSE;
+ }
+
+ if (!PropagateTo(xR, b)) {
+ if (AliLog::GetGlobalDebugLevel()>2) {
+ printf("Failed to propagate to X=%f for R=%f\n",xR,r);
+ Print("l");
+ }
+ return kFALSE;
+ }
+ double rcurr2 = xR*xR + GetY()*GetY();
+ if (TMath::Abs(rcurr2-rr)<kTiny || rr<kAlmost0) return kTRUE;
+ //
+ // two radii correspond to this X...
+ double pos[3]; GetXYZ(pos);
+ double phi = TMath::ATan2(pos[1],pos[0]);
+ if (!Rotate(phi)) {
+ if (AliLog::GetGlobalDebugLevel()>2) {
+ printf("Failed to rotate to %f to propagate to R=%f\n",phi,r);
+ Print("l");
+ }
+ return kFALSE;
+ }
+ if (++iter>8) {
+ if (AliLog::GetGlobalDebugLevel()>2) {
+ printf("Failed to propagate to R=%f after %d steps\n",r,iter);
+ Print("l");
+ }
+ return kFALSE;
+ }
+ }
+ return kTRUE;
+}
+
+
+//__________________________________________________________________________
+Bool_t KMCProbe::CorrectForMeanMaterial(const KMCLayer* lr, Bool_t inward)
+{
+ // printf("before at r=%.1f p=%.4f\n",lr->fR, P());
+ if (AliExternalTrackParam::CorrectForMeanMaterial(lr->fx2X0, inward ? lr->fXRho : -lr->fXRho, GetMass() , kTRUE)) {
+ // printf("after at r=%.1f p=%.4f\n",lr->fR, P());
+ return kTRUE;
+ }
+ AliDebug(2,Form("Failed to apply material correction, X/X0=%.4f", lr->fx2X0));
+ if (AliLog::GetGlobalDebugLevel()>1) Print();
+ return kFALSE;
+}
+
+/////////////////////////////////////////////////////////////////////////////
+ClassImp(KMCCluster)
+
+//_________________________________________________________________________
+KMCCluster::KMCCluster(KMCCluster &src)
+: TObject(src),
+ fY(src.fY),fZ(src.fZ),fX(src.fX),fPhi(src.fPhi)
+{}
+
+//__________________________________________________________________________
+KMCCluster& KMCCluster::operator=(const KMCCluster& src)
+{
+ if (this!=&src) {
+ TObject::operator=(src);
+ fY = src.fY;
+ fZ = src.fZ;
+ fX = src.fX;
+ fPhi = src.fPhi;
+ }
+ return *this;
+}
+
+//_________________________________________________________________________
+void KMCCluster::Print(Option_t *) const
+{
+ printf(" Local YZ = (%3.4lf,%3.4lf) | X=%3.4lf phi: %+.3f %s\n",fY,fZ,fX,fPhi,IsKilled()?"Killed":"");
+}
+
+/////////////////////////////////////////////////////////////////////////////
+ClassImp(KMCLayer)
+
+Double_t KMCLayer::fgDefEff = 1.0;
+//__________________________________________________________________________
+KMCLayer::KMCLayer(char *name) :
+ TNamed(name,name),fR(0),fx2X0(0),fPhiRes(0),fZRes(0),fEff(0),fIsDead(kFALSE),fType(-1),fActiveID(-1),fSig2EstD(999),fSig2EstZ(999),
+ fClCorr(),fClMC(),fClBg("KMCCluster",5), fTrCorr(), fTrMC("KMCProbe",5)
+{
+ Reset();
+}
+
+//__________________________________________________________________________
+void KMCLayer::Reset()
+{
+ fTrCorr.Reset();
+ fClCorr.Reset();
+ ResetMC();
+ fSig2EstD = fSig2EstZ = 999;
+ //
+}
+
+//__________________________________________________________________________
+KMCProbe* KMCLayer::AddMCTrack(KMCProbe* src)
+{
+ int ntr = GetNMCTracks();
+ KMCProbe* prb = 0;
+ if (src) prb = new(fTrMC[ntr]) KMCProbe(*src);
+ else prb = new(fTrMC[ntr]) KMCProbe();
+ if (!IsDead()) prb->ResetHit(GetActiveID());
+ return prb;
+}
+
+//__________________________________________________________________________
+void KMCLayer::Print(Option_t *opt) const
+{
+ printf("Lr%3d(A%3d) %10s R=%5.1f X2X0=%.3f XRho=%.3f SigY=%.4f SigZ=%.4f Eff:%4.2f\n",GetUniqueID(),fActiveID,GetName(), fR, fx2X0,fXRho,fPhiRes,fZRes,fEff);
+ TString opts = opt; opts.ToLower();
+ if (opts.Contains("c")) {
+ printf("Clusters: MC: %+7.4f:%+7.4f Ideal: %+7.4f:%+7.4f NBgCl: %3d NTrMC: %4d\n",fClMC.fY,fClMC.fZ, fClCorr.fY,fClCorr.fZ, GetNBgClusters(),GetNMCTracks());
+ }
+}
+
+/////////////////////////////////////////////////////////////////////////////
+Double_t KMCDetector::fgVtxConstraint[2]={-1,-1};
+
+ClassImp(KMCDetector)
+KMCDetector::KMCDetector() :
+TNamed("test_detector","detector"),
+ fNLayers(0),
+ fNActiveLayers(0),
+ fNActiveITSLayers(0),
+ fLastActiveLayer(-1),
+ fLastActiveITSLayer(-1),
+ fLastActiveLayerTracked(-1),
+ fBFieldG(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 ...
+ fdNdEtaCent(2300),
+ fMaxChi2Cl(25.),
+ fMaxNormChi2NDF(5.),
+ fMinITSHits(4),
+//
+ fMaxChi2ClSQ(4.), // precalulated internally
+ fMaxSeedToPropagate(300),
+//
+ fUseBackground(kFALSE),
+ fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
+ fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
+ fDensFactorEta(1.),
+//
+ fUpdCalls(0),
+ fHMCLrResidRPhi(0),
+ fHMCLrResidZ(0),
+ fHMCLrChi2(0),
+ fPattITS(0)
+{
+ //
+ // default constructor
+ //
+ // fLayers = new TObjArray();
+ RequireMaxChi2Cl(fMaxChi2Cl); // just to precalulate default square
+}
+
+KMCDetector::KMCDetector(char *name, char *title)
+ : TNamed(name,title),
+ fNLayers(0),
+ fNActiveLayers(0),
+ fNActiveITSLayers(0),
+ fLastActiveLayer(-1),
+ fLastActiveITSLayer(-1),
+ fLastActiveLayerTracked(-1),
+ fBFieldG(5.0),
+ 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 ...
+ fdNdEtaCent(2300),
+ fMaxChi2Cl(9.),
+ fMaxNormChi2NDF(5.),
+ fMinITSHits(4),
+ //
+ fMaxChi2ClSQ(3.), // precalulated internally
+ fMaxSeedToPropagate(50),
+ //
+ fUseBackground(kFALSE),
+ fBgYMin(1e6),fBgYMax(-1e6),fBgZMin(1e6),fBgZMax(-1e6),
+ fBgYMinTr(100),fBgYMaxTr(100),fBgZMinTr(100),fBgZMaxTr(100),fNBgLimits(0),
+ fDensFactorEta(1.),
+ //
+ fUpdCalls(0),
+ fHMCLrResidRPhi(0),
+ fHMCLrResidZ(0),
+ fHMCLrChi2(0),
+ fPattITS(0)
+{
+ //
+ // default constructor, that set the name and title
+ //
+ // fLayers = new TObjArray();
+}
+KMCDetector::~KMCDetector() { //
+ // virtual destructor
+ //
+ // delete fLayers;
+}
+
+void KMCDetector::InitMCWatchHistos()
+{
+ // init utility histos used for MC tuning
+ enum {kNBinRes=1000};
+ const double kMaxResidRPhi=1.0,kMaxResidZ=1.0,kMaxChi2=50;
+ int nlr = fNActiveITSLayers;
+ TString nam = "mc_residrhi";
+ TString tit = "MC $Delta Cl-Tr R#phi";
+ fHMCLrResidRPhi = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidRPhi,nlr,-0.5,nlr-0.5);
+ fHMCLrResidRPhi->GetXaxis()->SetTitle("cl-tr #Delta r#phi");
+ fHMCLrResidRPhi->Sumw2();
+ //
+ nam = "mc_residz";
+ tit = "MC $Delta Cl-Tr Z";
+ fHMCLrResidZ = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxResidZ,nlr,-0.5,nlr-0.5);
+ fHMCLrResidZ->GetXaxis()->SetTitle("cl-tr #Delta z");
+ fHMCLrResidZ->Sumw2();
+ //
+ nam = "mc_chi2";
+ tit = "MC #chi^{2} Cl-Tr Z";
+ fHMCLrChi2 = new TH2F(nam.Data(),tit.Data(),kNBinRes,-kMaxResidRPhi,kMaxChi2,nlr,-0.5,nlr-0.5);
+ fHMCLrChi2->GetXaxis()->SetTitle("cl-tr #chi^{2}");
+ fHMCLrChi2->Sumw2();
+ //
+ SetBit(kUtilHisto);
+}
+
+
+void KMCDetector::AddLayer(char *name, Float_t radius, Float_t x2X0, Float_t xrho, Float_t phiRes, Float_t zRes, Float_t eff) {
+ //
+ // Add additional layer to the list of layers (ordered by radius)
+ //
+
+ KMCLayer *newLayer = (KMCLayer*) fLayers.FindObject(name);
+
+ if (!newLayer) {
+ newLayer = new KMCLayer(name);
+ newLayer->fR = radius;
+ newLayer->fx2X0 = x2X0;
+ newLayer->fXRho = xrho;
+ newLayer->fPhiRes = phiRes;
+ newLayer->fZRes = zRes;
+ eff = TMath::Min(1.f,eff);
+ newLayer->fEff = eff <0 ? KMCLayer::GetDefEff() : eff;
+ newLayer->fActiveID = -2;
+ TString lname = name;
+ newLayer->fType = KMCLayer::kTypeNA;
+ if (lname.Contains("tpc")) newLayer->fType = KMCLayer::kTPC;
+ else if (lname.Contains("its")) newLayer->fType = KMCLayer::kITS;
+ if (lname.Contains("vertex")) newLayer->SetBit(KMCLayer::kBitVertex);
+ //
+ if (newLayer->fType==KMCLayer::kTypeNA) printf("Attention: the layer %s has undefined type\n",name);
+ //
+ newLayer->fIsDead = (newLayer->fPhiRes==RIDICULOUS && newLayer->fZRes==RIDICULOUS);
+ //
+ if (fLayers.GetEntries()==0)
+ fLayers.Add(newLayer);
+ else {
+ //
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ KMCLayer *l = (KMCLayer*)fLayers.At(i);
+ if (radius<l->fR) { fLayers.AddBefore(l,newLayer); break; }
+ if (radius>l->fR && (i+1)==fLayers.GetEntries() ) fLayers.Add(newLayer); // even bigger then last one
+ }
+ //
+ }
+ //
+ ClassifyLayers();
+ //
+ } else {
+ printf("Layer with the name %s does already exist\n",name);
+ }
+}
+
+//____________________________________________________________
+void KMCDetector::ClassifyLayers()
+{
+ // assign active Id's, etc
+ fLastActiveLayer = -1;
+ fLastActiveITSLayer = -1;
+ fNActiveLayers = 0;
+ fNActiveITSLayers = 0;
+ //
+ int nl = GetNLayers();
+ for (int il=0;il<nl;il++) {
+ KMCLayer* lr = GetLayer(il);
+ lr->SetUniqueID(il);
+ if (!lr->IsDead()) {
+ fLastActiveLayer = il;
+ lr->fActiveID = fNActiveLayers++;
+ if (lr->IsITS()) {
+ fLastActiveITSLayer = il;
+ fNActiveITSLayers++;
+ }
+ }
+ }
+ //
+ KMCProbe::SetNITSLayers(fNActiveITSLayers);
+}
+
+//____________________________________________________________
+void KMCDetector::KillLayer(char *name) {
+ //
+ // Marks layer as dead. Contribution only by Material Budget
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot mark as dead\n",name);
+ else {
+ tmp->fPhiRes = 999999;
+ tmp->fZRes = 999999;
+ tmp->fIsDead = kTRUE;
+ ClassifyLayers();
+ }
+}
+
+void KMCDetector::SetRadius(char *name, Float_t radius) {
+ //
+ // Set layer radius [cm]
+ //
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ //
+ if (!tmp) {
+ printf("Layer %s not found - cannot set radius\n",name);
+ } else {
+
+ Float_t tmpRadL = tmp->fx2X0;
+ Float_t tmpPhiRes = tmp->fPhiRes;
+ Float_t tmpZRes = tmp->fZRes;
+ Float_t tmpXRho = tmp->fXRho;
+ RemoveLayer(name); // so that the ordering is correct
+ AddLayer(name,radius,tmpRadL,tmpXRho,tmpPhiRes,tmpZRes);
+ }
+}
+
+Float_t KMCDetector::GetRadius(char *name) {
+ //
+ // Return layer radius [cm]
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get radius\n",name);
+ else
+ return tmp->fR;
+
+ return 0;
+}
+
+void KMCDetector::SetRadiationLength(char *name, Float_t x2X0) {
+ //
+ // Set layer material [cm]
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set layer material\n",name);
+ else {
+ tmp->fx2X0 = x2X0;
+ }
+}
+
+Float_t KMCDetector::GetRadiationLength(char *name) {
+ //
+ // Return layer radius [cm]
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get layer material\n",name);
+ else
+ return tmp->fx2X0;
+
+ return 0;
+
+}
+
+void KMCDetector::SetResolution(char *name, Float_t phiRes, Float_t zRes) {
+ //
+ // Set layer resolution in [cm]
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set resolution\n",name);
+ else {
+ tmp->fPhiRes = phiRes;
+ tmp->fZRes = zRes;
+ tmp->fIsDead = (zRes==RIDICULOUS && phiRes==RIDICULOUS);
+ ClassifyLayers();
+ }
+}
+
+Float_t KMCDetector::GetResolution(char *name, Int_t axis) {
+ //
+ // Return layer resolution in [cm]
+ // axis = 0: resolution in rphi
+ // axis = 1: resolution in z
+ //
+
+ KMCLayer *tmp = GetLayer(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get resolution\n",name);
+ else {
+ if (axis==0) return tmp->fPhiRes;
+ if (axis==1) return tmp->fZRes;
+ printf("error: axis must be either 0 or 1 (rphi or z axis)\n");
+ }
+ return 0;
+}
+
+void KMCDetector::SetLayerEfficiency(char *name, Float_t eff) {
+ //
+ // Set layer efficnecy (prop that his is missed within this layer)
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot set layer efficiency\n",name);
+ else {
+ tmp->fEff = eff;
+ }
+}
+
+Float_t KMCDetector::GetLayerEfficiency(char *name) {
+ //
+ // Get layer efficnecy (prop that his is missed within this layer)
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot get layer efficneicy\n",name);
+ else
+ return tmp->fEff;
+
+ return 0;
+
+}
+
+void KMCDetector::RemoveLayer(char *name) {
+ //
+ // Removes a layer from the list
+ //
+
+ KMCLayer *tmp = (KMCLayer*) fLayers.FindObject(name);
+ if (!tmp)
+ printf("Layer %s not found - cannot remove it\n",name);
+ else {
+ fLayers.Remove(tmp);
+ ClassifyLayers();
+ }
+}
+
+
+void KMCDetector::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]\n");
+
+ KMCLayer *tmp = 0;
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ tmp = (KMCLayer*)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->fR, tmp->fx2X0);
+ if (tmp->fPhiRes==RIDICULOUS)
+ printf(" - ");
+ else
+ printf("%3.0f ",tmp->fPhiRes*10000);
+ if (tmp->fZRes==RIDICULOUS)
+ printf(" -\n");
+ else
+ printf("%3.0f\n",tmp->fZRes*10000);
+ }
+}
+
+void KMCDetector::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);
+ KMCLayer *lastLayer = (KMCLayer*)fLayers.At(fLayers.GetEntries()-1); Double_t maxRad = lastLayer->fR;
+ gr->SetPointError(0,maxRad,maxRad);
+ gr->Draw("APE");
+
+
+ KMCLayer *tmp = 0;
+ for (Int_t i = fLayers.GetEntries()-1; i>=0; i--) {
+ tmp = (KMCLayer*)fLayers.At(i);
+
+
+ Double_t txtpos = tmp->fR;
+ 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->fR);
+ // 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 KMCDetector::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*)"IFCtpc", 77.8,0.01367, 0); // Inner Field cage (RS: set correct xrho for eloss)
+
+ // % Radiation Lengths ... Average per TPC row (i.e. total/159 )
+ Float_t x2X0PerRow = 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
+
+ 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,x2X0PerRow,0, phiResMean,zResMean);
+ else
+ AddLayer(Form("tpc_%d",k),rowRadius,x2X0PerRow,0); // non "active" row
+
+
+ }
+
+}
+
+void KMCDetector::RemoveTPC() {
+
+ // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-)
+ KMCLayer *tmp = 0;
+ for (Int_t i = 0; i<fLayers.GetEntries(); i++) {
+ tmp = (KMCLayer*)fLayers.At(i);
+ TString name(tmp->GetName());
+ if (name.Contains("tpc")) { RemoveLayer((char*)name.Data()); i--; }
+ }
+ RemoveLayer((char*)"IFC");
+
+}
+
+
+Double_t KMCDetector::ThetaMCS ( Double_t mass, Double_t x2X0, 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 ( x2X0 > 0 ) theta = 0.0136 * TMath::Sqrt(x2X0) / ( beta * momentum ) * (1+0.038*TMath::Log(x2X0)) ;
+ return (theta) ;
+}
+
+
+Double_t KMCDetector::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 KMCDetector::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 KMCDetector::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 KMCDetector::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 KMCDetector::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<0.01) return 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 KMCDetector::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) * fDensFactorEta ; // 2 eta ?
+ // note: surface of sphere is '4*pi*r^2'
+ // surface of cylinder is '2*pi*r* h'
+ return den ;
+}
+
+
+double KMCDetector::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.
+ double dens1 = OneEventHitDensity(multiplicity,radius);
+ if ( den < dens1 ) den = dens1;
+
+ return den ;
+}
+
+
+double KMCDetector::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 KMCDetector::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 KMCDetector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][20] ) 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(fBFieldG)) ;
+ Int_t kaonindex = (int)((ptK -0.1)*100.0 - 65.0*TMath::Abs(fBFieldG)) ;
+
+ if ( pionindex >= 20 ) 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 >= 20 ) 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 ;
+
+}
+
+KMCProbe* KMCDetector::PrepareKalmanTrack(double pt, double lambda, double mass, int charge, double phi, double x,double y, double z)
+{
+ // Prepare trackable Kalman track at the farthest position
+ //
+ // Set track parameters
+ // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis
+ fProbe.Reset();
+ fProbe.SetMass(mass);
+ KMCProbe* probe = new KMCProbe(fProbe);
+ double *trPars = (double*)probe->GetParameter();
+ double *trCov = (double*)probe->GetCovariance();
+ double xyz[3] = {x,y,z};
+ probe->Global2LocalPosition(xyz,phi);
+ probe->Set(xyz[0],phi,trPars,trCov);
+ trPars[KMCProbe::kY] = xyz[1];
+ trPars[KMCProbe::kZ] = xyz[2];
+ trPars[KMCProbe::kSnp] = 0; // track along X axis at the vertex
+ trPars[KMCProbe::kTgl] = TMath::Tan(lambda); // dip
+ trPars[KMCProbe::kPtI] = charge/pt; // q/pt
+ //
+ // put tiny errors to propagate to the outer-most radius
+ trCov[KMCProbe::kY2] = trCov[KMCProbe::kZ2] = trCov[KMCProbe::kSnp2] = trCov[KMCProbe::kTgl2] = trCov[KMCProbe::kPtI2] = 1e-20;
+ fProbe = *probe; // store original track
+ //
+ // propagate to last layer
+ fLastActiveLayerTracked = 0;
+ for (Int_t j=0; j<=fLastActiveLayer; j++) {
+ KMCLayer* lr = GetLayer(j);
+ lr->Reset();
+ //
+ if (!PropagateToLayer(probe,lr,1)) break;
+ if (!probe->CorrectForMeanMaterial(lr, kFALSE)) break;
+ //
+ lr->fClCorr.Set(probe->GetY(),probe->GetZ(), probe->GetX(), probe->GetAlpha());
+ if (!lr->IsDead()) fLastActiveLayerTracked = j;
+ }
+ probe->ResetCovMat();// reset cov.matrix
+ printf("Last active layer trracked: %d (out of %d)\n",fLastActiveLayerTracked,fLastActiveLayer);
+ //
+ return probe;
+}
+
+
+
+TGraph * KMCDetector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) {
+ //
+ // returns the momentum resolution
+ //
+
+ TGraph *graph = new TGraph(20, 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 * KMCDetector::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 ( 20, fTransMomenta, fResolutionRPhi ) ;
+ graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ;
+ graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)") ;
+ } else {
+ graph = new TGraph ( 20, fTransMomenta, fResolutionZ ) ;
+ 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->SetLineWidth(linewidth);
+ graph->SetLineColor(color);
+ graph->SetMarkerColor(color);
+
+ return graph;
+
+}
+
+
+TGraph * KMCDetector::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[20];
+
+ 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++) {
+ KMCLayer *l = (KMCLayer*)fLayers.At(i);
+ if (!l->IsDead() && l->fR>0) {
+ layerRadius[count] = l->fR;
+ layerThickness[count] = l->fx2X0;
+ if (axis==0) {
+ layerResolution[count] = l->fPhiRes;
+ } else {
+ layerResolution[count] = l->fZRes;
+ }
+ 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 < 20 ; 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 ( 20, 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 * KMCDetector::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[20]; // with chosen particle mass
+ Double_t kaonEfficiency[20], pionEfficiency[20], d0efficiency[20];
+ Double_t partEfficiency[2][20];
+
+ 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 < 20 ; 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 < 20 ; j++ ) {
+ d0efficiency[j] = D0IntegratedEfficiency(fTransMomenta[j],partEfficiency);
+ }
+ } else {
+ for ( Int_t j = 0 ; j < 20 ; j++ ) {
+ particleEfficiency[j] = fEfficiency[2][j]* AcceptanceOfTpcAndSi;
+ // NOTE: Decay factor (see kaon) should be included to be realiable
+ }
+ }
+
+ for ( Int_t j = 0 ; j < 20 ; j++ ) {
+ pionEfficiency[j] *= 100;
+ kaonEfficiency[j] *= 100;
+ d0efficiency[j] *= 100;
+ particleEfficiency[j] *= 100;
+ }
+
+ TGraph * graph = 0;
+ if (particle==0) {
+ graph = new TGraph ( 20, fTransMomenta, particleEfficiency ) ; // choosen mass
+ graph->SetLineWidth(1);
+ } else if (particle==1) {
+ graph = new TGraph ( 20, fTransMomenta, pionEfficiency ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==2) {
+ graph = new TGraph ( 20, fTransMomenta, kaonEfficiency ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==3) {
+ graph = new TGraph ( 20, 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 * KMCDetector::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[20]; // with chosen particle mass
+ Double_t kaonFake[20], pionFake[20];
+ Double_t partFake[2][20];
+
+ 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 < 20 ; 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 < 20 ; j++ ) {
+ particleFake[j] = fFake[2][j];
+ // NOTE: Decay factor (see kaon) should be included to be realiable
+ }
+ }
+
+ for ( Int_t j = 0 ; j < 20 ; j++ ) {
+ pionFake[j] *= 100;
+ kaonFake[j] *= 100;
+ particleFake[j] *= 100;
+ }
+
+ TGraph * graph = 0;
+ if (particle==0) {
+ graph = new TGraph ( 20, fTransMomenta, particleFake ) ; // choosen mass
+ graph->SetLineWidth(1);
+ } else if (particle==1) {
+ graph = new TGraph ( 20, fTransMomenta, pionFake ) ;
+ graph->SetLineWidth(1);
+ } else if (particle ==2) {
+ graph = new TGraph ( 20, 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* KMCDetector::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* KMCDetector::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 KMCDetector::MakeAliceAllNew(Bool_t flagTPC,Bool_t flagMon, int setVer) {
+
+ // All New configuration with X0 = 0.3 and resolution = 4 microns
+
+ AddLayer((char*)"bpipe_its",2.0,0.0022, 0.092); // beam pipe, 0.5 mm Be
+ AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation
+ if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
+ printf("vertex will work as constraint: %.4f %.4f\n",fgVtxConstraint[0],fgVtxConstraint[1]);
+ SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
+ }
+ //
+ // new ideal Pixel properties?
+ Double_t x0 = 0.0050;
+ Double_t resRPhi = 0.0006;
+ Double_t resZ = 0.0006;
+ Double_t xrho = 0.0116; // assume 0.5mm of silicon for eloss
+ //
+ if (flagMon) {
+ x0 *= 3./5.;
+ xrho *=3./5.;
+ resRPhi = 0.0004;
+ resZ = 0.0004;
+ }
+
+ double sclD = 1;
+ double sclZ = 1;
+
+ // default all new
+ if (setVer<=0) {
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
+ }
+ else if (setVer==1) {
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 2.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 3.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its5", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ //
+ /*
+ UInt_t patt[] = {
+ KMCTrackSummary::Bits(1,1,1),
+ KMCTrackSummary::Bits(0,0,0,1,1),
+ KMCTrackSummary::Bits(0,0,0,0,0,1,1)
+ };
+ RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
+ */
+ }
+ else if (setVer==2) {
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 20.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its4", 22.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its5", 33.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its6", 43.0 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ AddLayer((char*)"its7", 43.6 , x0, xrho, resRPhi*sclD, resZ*sclZ);
+ //
+ /*
+ UInt_t patt[] = {
+ KMCTrackSummary::Bits(1,1,1,1),
+ KMCTrackSummary::Bits(0,0,0,0,1,1),
+ KMCTrackSummary::Bits(0,0,0,0,0,0,1,1,1)
+ };
+ RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
+ */
+ }
+ /*
+ // last 2 layers strips
+ double resRPStr=0.0020, resZStr = 0.0830;
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6", 39.6 , x0, xrho, resRPStr, resZStr);
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPStr, resZStr);
+ */
+ //*/
+ /*
+ // resolution scaled with R as res(2.2) * R/2.2
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
+ KMCLayer* lr0 = 0;
+ for (int i=0;i<GetNActiveITSLayers();i++) {
+ KMCLayer *lr = GetActiveLayer(i);
+ if (lr->IsVertex()) continue;
+ if (lr0==0) {
+ lr0=lr;
+ // continue;
+ }
+ double scl = 5*lr->GetRadius()/lr0->GetRadius();
+ SetResolution((char*)lr->GetName(), resRPhi*scl, resZ*scl*4);
+ }
+ */
+ /*
+ // 1st 2 layers double
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its1a", 2.8 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2a", 4.2 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its3a", 6.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
+ */
+ /*
+ // last 4 layers doubled
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4a", 12.8 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its5", 23.5 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5a", 23.9 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its6", 39.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6a", 40.0 , x0, xrho, resRPhi, resZ);
+
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its7a", 43.4 , x0, xrho, resRPhi, resZ);
+ */
+
+ /* //last 3 lr together
+ AddLayer((char*)"its1", 2.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its2", 3.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its3", 6.8 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its4", 12.4 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its5", 42.2 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its6", 42.6 , x0, xrho, resRPhi, resZ);
+ AddLayer((char*)"its7", 43.0 , x0, xrho, resRPhi, resZ);
+ */
+ if (flagTPC) {
+ AddTPC(0.1,0.1); // TPC
+ }
+ PrintLayout();
+}
+
+void KMCDetector::MakeAliceCurrent(Bool_t flagTPC, Int_t AlignResiduals) {
+
+ // 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
+
+ // 1.48e-01, 2.48e-01,2.57e-01, 1.34e-01, 3.34e-01,3.50e-01, 2.22e-01, 2.38e-01,2.25e-01};
+
+ AddLayer((char*)"vertex_its", 0, 0); // dummy vertex for matrix calculation
+ if (fgVtxConstraint[0]>0 && fgVtxConstraint[1]>0) {
+ printf("vertex will work as constraint: %.4f %.4f",fgVtxConstraint[0],fgVtxConstraint[1]);
+ SetResolution((char*)"vertex_its",fgVtxConstraint[0],fgVtxConstraint[1]);
+ }
+ //
+ AddLayer((char*)"bpipe_its",2.94,0.0022, 1.48e-01); // beam pipe
+ AddLayer((char*)"tshld1_its",11.5,0.0065, 1.34e-01); // Thermal shield // 1.3% /2
+ AddLayer((char*)"tshld2_its",31.0,0.0108, 2.22e-01); // Thermal shield // 1.3% /2
+ //
+ if (flagTPC) {
+ AddTPC(0.1,0.1); // TPC
+ }
+ // Adding the ITS - current configuration
+
+ if (AlignResiduals==0) {
+ AddLayer((char*)"spd1_its", 3.9, 0.0114, 2.48e-01, 0.0012, 0.0130);
+ AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, 0.0012, 0.0130);
+ AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, 0.0035, 0.0025);
+ AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, 0.0035, 0.0025);
+ AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, 0.0020, 0.0830);
+ AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 0.0020, 0.0830);
+ /*
+ UInt_t patt[] = {
+ KMCTrackSummary::Bits(1,1),
+ KMCTrackSummary::Bits(0,0,1,1),
+ KMCTrackSummary::Bits(0,0,0,0,1,1)
+ };
+ RequirePattern( patt, sizeof(patt)/sizeof(UInt_t));
+ */
+ } 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_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030),
+ TMath::Sqrt(0.0130*0.0130+0.0050*0.0050));
+ AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0100*0.0100),
+ TMath::Sqrt(0.0025*0.0025+0.0050*0.0050));
+ AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020),
+ TMath::Sqrt(0.0830*0.0830+0.1000*0.1000));
+ AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, 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_its", 3.9, 0.0114, 2.48e-01, 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_ots", 7.6, 0.0114, 2.57e-01, 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_ots",15.0, 0.0113, 3.34e-01, 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_its",23.9, 0.0126, 3.50e-01, 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_its",38.0, 0.0083, 2.38e-01, 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_its",43.0, 0.0086, 2.25e-01, 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_its", 3.9, 0.0114, 2.48e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.000*0.000));
+ AddLayer((char*)"spd2_its", 7.6, 0.0114, 2.57e-01, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008),
+ TMath::Sqrt(0.0130*0.0130+0.000*0.000));
+ AddLayer((char*)"sdd1_its",15.0, 0.0113, 3.34e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.000*0.000));
+ AddLayer((char*)"sdd2_its",23.9, 0.0126, 3.50e-01, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010),
+ TMath::Sqrt(0.0025*0.0025+0.000*0.000));
+ AddLayer((char*)"ssd1_its",38.0, 0.0083, 2.38e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.000*0.000));
+ AddLayer((char*)"ssd2_its",43.0, 0.0086, 2.25e-01, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010),
+ TMath::Sqrt(0.0830*0.0830+0.000*0.000));
+ }
+
+}
+
+
+void KMCDetector::MakeStandardPlots(Bool_t add, Int_t color, Int_t linewidth,Bool_t onlyPionEff) {
+ //
+ // Produces the standard performace plots
+ //
+
+ 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();
+ TGraph *eff = GetGraphRecoEfficiency(1,color,linewidth);
+ eff->SetTitle("Efficiencies");
+ eff->Draw("AL");
+ if (!onlyPionEff) {
+ GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
+ GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
+ }
+ c1->cd(2); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogy(); gPad->SetLogx();
+ GetGraphMomentumResolution(color,linewidth)->Draw("AL");
+
+ c1->cd(3); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogx();
+ GetGraphPointingResolution(0,color,linewidth)->Draw("AL");
+
+ c1->cd(4); gPad->SetGridx(); gPad->SetGridy();
+ gPad->SetLogx();
+ GetGraphPointingResolution(1,color,linewidth)->Draw("AL");
+
+ } else {
+
+ TVirtualPad *c1 = gPad->GetMother();
+
+ c1->cd(1);
+ GetGraphRecoEfficiency(1,color,linewidth)->Draw("L");
+ if (!onlyPionEff) {
+ GetGraphRecoEfficiency(2,color,linewidth)->Draw("L");
+ GetGraphRecoEfficiency(3,color,linewidth)->Draw("L");
+ }
+ c1->cd(2); GetGraphMomentumResolution(color,linewidth)->Draw("L");
+
+ c1->cd(3); GetGraphPointingResolution(0,color,linewidth)->Draw("L");
+
+ c1->cd(4); GetGraphPointingResolution(1,color,linewidth)->Draw("L");
+
+ }
+
+}
+
+void KMCDetector::ApplyMS(KMCProbe* trc, double x2X0) const
+{
+ // simulate random modification of track params due to the MS
+ if (x2X0<=0) return;
+ double alpha = trc->GetAlpha(); // store original alpha
+ double mass = trc->GetMass();
+ //
+ double snp = trc->GetSnp();
+ double dip = trc->GetTgl();
+ Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));
+ x2X0 *= angle;
+ //
+ static double covCorr[15],covDum[21]={0};
+ static double mom[3],pos[3];
+ double *cov = (double*) trc->GetCovariance();
+ memcpy(covCorr,cov,15*sizeof(double));
+ trc->GetXYZ(pos);
+ trc->GetPxPyPz(mom);
+ double pt2 = mom[0]*mom[0]+mom[1]*mom[1];
+ double pt = TMath::Sqrt(pt2);
+ double ptot2 = pt2 + mom[2]*mom[2];
+ double ptot = TMath::Sqrt(ptot2);
+ double beta = ptot/TMath::Sqrt(ptot2 + mass*mass);
+ double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);
+ //
+ // a la geant
+ double phiSC = gRandom->Rndm()*TMath::Pi();
+ double thtSC = gRandom->Gaus(0,1.4142*sigth);
+ // printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);
+ double sn = TMath::Sin(thtSC);
+ double dx = sn*TMath::Sin(phiSC);
+ double dy = sn*TMath::Cos(phiSC);
+ double dz = TMath::Cos(thtSC);
+ double v[3];
+ // printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);
+ for (int i=3;i--;) mom[i] /= ptot;
+ double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+ if (!IsZero(pt)) {
+ double pd1 = mom[0]/vmm;
+ double pd2 = mom[1]/vmm;
+ v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;
+ v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;
+ v[2] = -vmm*dx + mom[2]*dz;
+ }
+ else {
+ v[0] = dx;
+ v[1] = dy;
+ v[2] = dz*TMath::Sign(1.,mom[2]);
+ }
+ double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+ // printf("before :%+e %+e %+e || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2], sigth, x2X0, pt, beta);
+ // trc->Print();
+ // direction cosines -> p
+ for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;
+ // printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);
+ trc->Set(pos,mom,covDum,trc->Charge());
+ //
+ trc->Rotate(alpha);
+ memcpy(cov,covCorr,15*sizeof(double));
+ //
+}
+
+//________________________________________________________________________________
+Bool_t KMCDetector::TransportKalmanTrackWithMS(KMCProbe *probTr, int maxLr) const
+{
+ // Transport track till layer maxLr, applying random MS
+ //
+ for (Int_t j=0; j<maxLr; j++) {
+ KMCLayer* lr0 = (KMCLayer*)fLayers.At(j);
+ KMCLayer* lr = (KMCLayer*)fLayers.At(j+1);
+ if (!lr0->IsVertex()) {
+ ApplyMS(probTr,lr0->fx2X0); // apply MS
+ if (!probTr->CorrectForMeanMaterial(lr0,kFALSE)) return kFALSE;
+ }
+ //
+ if (!PropagateToLayer(probTr,lr,1)) return kFALSE;
+ // store randomized cluster local coordinates and phi
+ lr->ResetBgClusters();
+ double rz,ry;
+ gRandom->Rannor(rz,ry);
+ lr->GetMCCluster()->Set(probTr->GetY()+ry*lr->GetPhiRes(),probTr->GetZ()+rz*lr->GetZRes(),
+ probTr->GetX(), probTr->GetAlpha() );
+ //
+ }
+ //
+ return kTRUE;
+}
+
+//________________________________________________________________________________
+Bool_t KMCDetector::SolveSingleTrack(Double_t mass, Double_t pt, Double_t eta, TObjArray* sumArr,int nMC, int offset)
+{
+ // analityc and fullMC (nMC trials) evaluaion of tracks with given kinematics.
+ // the results are filled in KMCTrackSummary objects provided via summArr array
+ //
+ int progressP = 10;//0; // print progress in percents
+ //
+ progressP = int(nMC*0.01*progressP);
+
+ if (!SolveSingleTrackViaKalman(mass,pt,eta)) return kFALSE;
+ //
+ // Store non-updated track errors of inward propagated seed >>>>>>>>
+ int maxLr = fLastActiveITSLayer + offset;
+ if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
+ KMCProbe probeTmp = fProbe; // original probe at vertex
+ KMCLayer* lr = 0;
+ for (Int_t j=1; j<=maxLr; j++) {
+ lr = GetLayer(j);
+ // printf("Here0: %d\n",j);
+ if (!PropagateToLayer(&probeTmp,lr,1)) return 0;
+ if (j!=maxLr) if (!probeTmp.CorrectForMeanMaterial(lr, kFALSE)) return 0;
+ // printf("Prelim. Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(probeTmp.GetSigmaY2()),TMath::Sqrt(probeTmp.GetSigmaZ2()));
+ }
+ for (Int_t j=maxLr; j>0; j--) {
+ lr = GetLayer(j);
+ // printf("Here1: %d\n",j);
+ if (j!=maxLr) if (!PropagateToLayer(&probeTmp,lr,-1)) return 0;
+ lr->fSig2EstD = probeTmp.GetSigmaY2();
+ lr->fSig2EstZ = probeTmp.GetSigmaZ2();
+ // probeTmp.Print("l");
+ printf("Natural Err at lr:%8s | %7.3f %7.3f\n",lr->GetName(),TMath::Sqrt(lr->fSig2EstD),TMath::Sqrt(lr->fSig2EstZ));
+ if (!probeTmp.CorrectForMeanMaterial(lr, kTRUE)) return 0;
+ }
+ // Store non-updated track errors of inward propagated seed <<<<<<<<
+ //
+ int nsm = sumArr ? sumArr->GetEntriesFast() : 0;
+ KMCLayer* vtx = GetLayer(0);
+ //
+ for (int i=0;i<nsm;i++) {
+ KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(i);
+ if (!tsm) continue;
+ tsm->SetRefProbe( GetProbeTrack() ); // attach reference track (generated)
+ tsm->SetAnProbe( vtx->GetAnProbe() ); // attach analitycal solution
+ }
+ //
+ TStopwatch sw;
+ sw.Start();
+ for (int it=0;it<nMC;it++) {
+ printf("ev: %d\n",it);
+ SolveSingleTrackViaKalmanMC(offset);
+ KMCProbe* trc = vtx->GetWinnerMCTrack();
+ vtx->GetMCTracks()->Print();
+ if (progressP==1 || (progressP>0 && (it%progressP)==0)) {
+ printf("%d%% done |",it*100/nMC);
+ sw.Stop(); sw.Print(); sw.Start(kFALSE);
+ }
+ for (int ism=nsm;ism--;) { // account the track in each of summaries
+ KMCTrackSummary* tsm = (KMCTrackSummary*)sumArr->At(ism);
+ if (!tsm) continue;
+ tsm->AddUpdCalls(GetUpdCalls());
+ tsm->AddTrack(trc);
+ }
+ }
+ //
+ sw.Stop();
+ printf("Total time: "); sw.Print();
+ return kTRUE;
+}
+
+//________________________________________________________________________________
+Int_t KMCDetector::GetLayerID(int actID) const
+{
+ // find physical layer id from active id
+ if (actID<0 || actID>fNActiveLayers) return -1;
+ int start = actID<fNActiveITSLayers ? fLastActiveITSLayer : fLastActiveLayer;
+ for (int i=start+1;i--;) {
+ if (GetLayer(i)->GetActiveID()==actID) return i;
+ }
+ return -1;
+}
+
+//________________________________________________________________________________
+KMCProbe* KMCDetector::KalmanSmooth(int actLr, int actMin,int actMax) const
+{
+ // estimate kalman smoothed track params at given active lr
+ // from fit at layers actMin:actMax (excluding actLr)
+ // SolveSingleTrackViaKalman must have been called before
+ //
+ if (actMin>actMax) swap(actMin,actMax);
+ if (actMax>=fNActiveLayers) actMax = fNActiveLayers-1;
+ int nlrfit = actMax-actMin;
+ if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
+ if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
+ static KMCProbe iwd,owd;
+ //
+ // find phisical layer id's
+ int pLr = GetLayerID(actLr);
+ int pMin = GetLayerID(actMin);
+ int pMax = GetLayerID(actMax);
+ //
+ // printf(">>> %d %d %d\n",pLr, pMin,pMax);
+ Bool_t useIwd=kFALSE, useOwd=kFALSE;
+ if (pLr<pMax) { // need inward piece
+ iwd = GetLayer(pMax)->fTrCorr;
+ iwd.ResetCovMat();
+ iwd.GetHitsPatt() = 0;
+ for (int i=pMax;i>=pLr;i--) {
+ KMCLayer* lr = GetLayer(i);
+ // printf("IWD %d\n",i);
+ if (!lr->IsDead() && i!=pLr && i>=pMin) if (!UpdateTrack(&iwd,lr,&lr->fClCorr)) return 0;
+ if (i!=pLr) {
+ if (!iwd.CorrectForMeanMaterial(lr,kTRUE)) return 0; // correct for materials of this layer
+ if (!PropagateToLayer(&iwd,GetLayer(i-1),-1)) return 0; // propagate to next layer
+ }
+ // printf("IWD%d: ",i); iwd.Print("l");
+ }
+ useIwd = kTRUE;
+ }
+ if (pLr>pMin) { // need outward piece
+ owd = GetLayer(pMin)->fTrCorr;
+ owd.ResetCovMat();
+ owd.GetHitsPatt() = 0;
+ for (int i=pMin;i<=pLr;i++) {
+ KMCLayer* lr = GetLayer(i);
+ // printf("OWD %d\n",i);
+ if (!lr->IsDead() && i!=pLr && i<=pMax) if (!UpdateTrack(&owd,lr,&lr->fClCorr)) return 0;
+ if (i!=pLr) {
+ if (!owd.CorrectForMeanMaterial(lr,0)) return 0; // correct for materials of this layer
+ if (!PropagateToLayer(&owd,GetLayer(i+1), 1)) return 0; // propagate to next layer
+ }
+ // printf("OWD%d: ",i); owd.Print("l");
+ }
+ useOwd = kTRUE;
+ }
+ //
+ // was this extrapolation outside the fit range?
+ if (!useIwd) return (KMCProbe*)&owd;
+ if (!useOwd) return (KMCProbe*)&iwd;
+ //
+ // weight both tracks
+ if (!iwd.Propagate(owd.GetAlpha(),owd.GetX(),fBFieldG)) return 0;
+ double meas[2] = {owd.GetY(),owd.GetZ()};
+ double measErr2[3] = {owd.GetSigmaY2(), owd.GetSigmaZY(), owd.GetSigmaZ2()};
+ // printf("Weighting\n");
+ // owd.Print("l");
+ // iwd.Print("l");
+ if (!iwd.Update(meas,measErr2)) return 0;
+ iwd.GetHitsPatt() |= owd.GetHitsPatt();
+
+ // printf("->\n");
+ // iwd.Print("l");
+
+ return (KMCProbe*)&iwd;
+ //
+}
+
+//________________________________________________________________________________
+KMCProbe* KMCDetector::KalmanSmoothFull(int actLr, int actMin,int actMax) const
+{
+ // estimate kalman smoothed track params at given active lr
+ // from fit at layers actMin:actMax (excluding actLr)
+ // SolveSingleTrackViaKalman must have been called before
+ //
+ static TClonesArray prediction("KMCProbe",10);
+ static TClonesArray update("KMCProbe",10);
+ static KMCProbe res;
+ //
+ if (actMin>actMax) swap(actMin,actMax);
+ int nlrfit = actMax-actMin;
+ if (actLr>=actMin && actLr<=actMax) nlrfit-=1;
+ if (nlrfit<2) {AliInfo("Need a least 2 active layers in the fit"); return 0;}
+ //
+ // find phisical layer id's
+ int pLr = GetLayerID(actLr);
+ int pMin = GetLayerID(actMin);
+ int pMax = GetLayerID(actMax);
+ //
+ int dir=0,dirInt=0;
+ if (pLr<=pMin) dir=-1; // inward extrapolation
+ else if (pLr>=pMax) dir= 1; // outward extrapolation
+ else if (actMax-actLr >= actLr-actMin) dirInt = -1; // inward interpolation (the test point is closer to inner layer)
+ else dirInt = 1; // outward interpolation (the test point is closer to outer layer)
+ //
+ if (dir!=0) { // no sens to do smoothing: simple Kalman filtering extrapolation
+ int start = dir<0 ? pMax : pMin;
+ res = GetLayer(start)->fTrCorr;
+ res.ResetCovMat();
+ KMCLayer* lr = 0;
+ for (int i=(dir<0?pMax:pMin); i!=pLr; i+=dir) { // track till nearest layer to pLr
+ lr = GetLayer(i);
+ if (!lr->IsDead() && !(i<pMin ||i>pMax)) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range
+ if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer
+ if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer
+ }
+ if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this nearest layer
+ if (!PropagateToLayer(&res,GetLayer(pLr), dir)) return 0; // propagate to test layer
+ return (KMCProbe*)&res;
+ }
+ //
+ // too bad, need to do real filtering
+ //
+ int start = dirInt<0 ? pMax : pMin;
+ int stop = dirInt<0 ? pMin-1 : pMax+1;
+ res = GetLayer(start)->fTrCorr;
+ res.ResetCovMat();
+ KMCLayer* lr = 0;
+ int count = 0;
+ for (int i=start; i!=stop; i+=dirInt) { // track in full range, storing updates and predictions
+ new(prediction[count]) KMCProbe(res);
+ lr = GetLayer(i);
+ if (!lr->IsDead() && i!=pLr) if (!UpdateTrack(&res,lr,&lr->fClCorr)) return 0; // update only with layers in fit range
+ new(update[count]) KMCProbe(res);
+ if (!res.CorrectForMeanMaterial(lr,dir<0 ? kTRUE:kFALSE)) return 0; // correct for materials of this layer
+ if (!PropagateToLayer(&res,GetLayer(i+dir),dir)) return 0; // propagate to next layer
+ count++;
+ }
+ return (KMCProbe*)&res;
+ //
+}
+
+//________________________________________________________________________________
+Bool_t KMCDetector::SolveSingleTrackViaKalman(Double_t mass, Double_t pt, Double_t eta)
+{
+ // analytical estimate of tracking resolutions
+ // fProbe.SetUseLogTermMS(kTRUE);
+ //
+ if (fMinITSHits>fNActiveITSLayers) {fMinITSHits = fNActiveITSLayers; printf("Redefined request of min N ITS hits to %d\n",fMinITSHits);}
+ if (TMath::Abs(eta)<1e-3) fDensFactorEta = 1.;
+ else {
+ fDensFactorEta = TMath::Tan( 2.*TMath::ATan(TMath::Exp(-TMath::Abs(eta))) );
+ fDensFactorEta = 1./TMath::Sqrt( 1. + 1./fDensFactorEta/fDensFactorEta);
+ }
+ double lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-eta));
+ KMCProbe* probe = PrepareKalmanTrack(pt,lambda,mass,-1);
+ if (!probe) return kFALSE;
+ //
+ KMCLayer *lr = 0;
+ //
+ //
+ // Start the track fitting --------------------------------------------------------
+ //
+ // Back-propagate the covariance matrix along the track.
+ // Kalman loop over the layers
+ //
+ KMCProbe* currTr = 0;
+ lr = (KMCLayer*)fLayers.At(fLastActiveLayerTracked);
+ lr->fTrCorr = *probe;
+ delete probe; // rethink...
+ //
+ for (Int_t j=fLastActiveLayerTracked; j--; ) { // Layer loop
+ //
+ KMCLayer *lrP = lr;
+ lr = (KMCLayer*)fLayers.At(j);
+ //
+ lr->fTrCorr = lrP->fTrCorr;
+ currTr = &lr->fTrCorr;
+ currTr->ResetHit(lrP->GetActiveID());
+ //
+ // if there was a measurement on prev layer, update the track
+ if (!lrP->IsDead()) { // include measurement
+ KMCCluster cl(currTr->GetY(),currTr->GetZ(), currTr->GetX(), currTr->GetAlpha());
+ if (!UpdateTrack(currTr,lrP,&cl)) return kFALSE;
+ }
+ if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) return kFALSE; // correct for materials of this layer
+ if (!PropagateToLayer(currTr,lr,-1)) return kFALSE; // propagate to current layer
+ //
+ } // end loop over layers
+ //
+ return kTRUE;
+}
+
+//____________________________________________________________
+Bool_t KMCDetector::SolveSingleTrackViaKalmanMC(int offset)
+{
+ // MC estimate of tracking resolutions/effiencies. Requires that the SolveSingleTrackViaKalman
+ // was called before, since it uses data filled by this method
+ //
+ // The MC tracking will be done starting from fLastActiveITSLayer + offset (before analytical estimate will be used)
+ //
+ // At this point, the fProbe contains the track params generated at vertex.
+ // Clone it and propagate to target layer to generate hit positions affected by MS
+ //
+ fUpdCalls = 0.;
+ KMCProbe *currTrP=0,*currTr=0;
+ int maxLr = fLastActiveITSLayer + offset;
+ if (maxLr >= fLastActiveLayerTracked-1) maxLr = fLastActiveLayerTracked;
+ ResetMCTracks(maxLr);
+ KMCLayer* lr = (KMCLayer*)fLayers.At(maxLr);
+ currTr = lr->AddMCTrack(&fProbe); // start with original track at vertex
+ //
+ if (!TransportKalmanTrackWithMS(currTr, maxLr)) return kFALSE; // transport it to outermost layer where full MC is done
+ //
+ if (fLastActiveITSLayer<fLastActiveLayerTracked) { // prolongation from TPC
+ // start from correct track propagated from above till maxLr
+ double *covMS = (double*)currTr->GetCovariance();
+ const double *covIdeal =lr->fTrCorr.GetCovariance();
+ for (int i=15;i--;) covMS[i] = covIdeal[i];
+ }
+ else { // ITS SA: randomize the starting point
+ // double *pars = (double*)currTr->GetParameter();
+ // pars[0] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaY2()));
+ // pars[1] += gRandom->Gaus(0,TMath::Sqrt(currTr->GetSigmaZ2()));
+ //
+ currTr->ResetCovMat();
+ /*
+ double *trCov = (double*)currTr->GetCovariance();
+ double *trPars = (double*)currTr->GetParameter();
+ const double kLargeErr2PtI = 0.3*0.3;
+ trCov[14] = TMath::Max(trCov[14],kLargeErr2PtI*trPars[4]*trPars[4]);
+ */
+ }
+ //
+ for (Int_t j=maxLr; j--; ) { // Layer loop
+ //
+ KMCLayer *lrP = lr;
+ lr = (KMCLayer*)fLayers.At(j);
+ int ntPrev = lrP->GetNMCTracks();
+ //
+ if (lrP->IsDead()) { // for passive layer just propagate the copy of all tracks of prev layer >>>
+ for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
+ currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
+ currTr = lr->AddMCTrack( currTrP );
+ if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
+ if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer
+ }
+ continue;
+ } // treatment of dead layer <<<
+ //
+ if (lrP->IsTPC()) { // we don't consider bg hits in TPC, just update with MC cluster
+ for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
+ currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
+ currTr = lr->AddMCTrack( currTrP );
+ if (!UpdateTrack(currTr, lrP, lrP->GetMCCluster(), kTRUE)) {currTr->Kill(); continue;} // update with correct MC cl.
+ if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill(); continue;} // correct for materials of prev. layer
+ if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill(); continue;} // propagate to current layer
+ }
+ continue;
+ } // treatment of ideal (TPC?) layer <<<
+ //
+ // active layer under eff. study (ITS?): propagate copy of every track to MC cluster frame (to have them all in the same frame)
+ // and calculate the limits of bg generation
+ KMCCluster* clMC = lrP->GetMCCluster();
+ if (lrP->GetLayerEff()<gRandom->Rndm()) clMC->Kill(); // simulate inefficiency
+ ResetSearchLimits();
+ int nseeds = 0;
+ for (int itrP=ntPrev;itrP--;) { // loop over all tracks from previous layer
+ currTrP = lrP->GetMCTrack(itrP); if (currTrP->IsKilled()) continue;
+ currTr = lr->AddMCTrack( currTrP );
+ if (!currTr->PropagateToCluster(clMC,fBFieldG)) {currTr->Kill(); continue;} // go to MC cluster
+ if ( !(currTr->GetNITSHits()>0 && currTr->GetNITSHits()==currTr->GetNFakeITSHits()) ) UpdateSearchLimits(currTr, lrP); // RS
+ nseeds++;
+ }
+ //
+ // printf("%3d seeds\n",nseeds);
+ if (fUseBackground && lrP->IsITS()) GenBgClusters(lrP); // generate background hits
+ //
+ ntPrev = lr->GetNMCTracks();
+ for (int itr=ntPrev;itr--;) { // loop over all tracks PROPAGATED from previous layer to clusters frame on previous layer
+ currTrP = lr->GetMCTrack(itr); // this is a seed from prev layer. The new clusters are attached to its copies, the seed itself
+ // will be propagated w/o cluster update if it does not violate requested "reconstructed" track settings
+ if (currTrP->IsKilled()) continue;
+ //printf("Check %d %p %d\n",itr,currTrP,currTrP->GetUniqueID()); currTrP->Print();
+ CheckTrackProlongations(currTrP, lr, lrP);
+ if (NeedToKill(currTrP)) currTrP->Kill(); // kill track which was not updated at lrP
+ //currTrP->Kill(); // kill track which was not updated at lrP
+ }
+ //
+ lr->GetMCTracks()->Sort();
+ int ntTot = lr->GetNMCTracks(); // propagate max amount of allowed tracks to current layer
+ if (ntTot>fMaxSeedToPropagate && fMaxSeedToPropagate>0) {
+ for (int itr=ntTot;itr>=fMaxSeedToPropagate;itr--) lr->GetMCTracks()->RemoveAt(itr);
+ ntTot = fMaxSeedToPropagate;
+ }
+ //
+ for (int itr=ntTot;itr--;) {
+ currTr = lr->GetMCTrack(itr);
+ if (!currTr->CorrectForMeanMaterial(lrP,kTRUE)) {currTr->Kill();continue;} // correct for materials of prev. layer
+ if (!PropagateToLayer(currTr,lr,-1)) {currTr->Kill();continue;} // propagate to current layer
+ }
+ AliDebug(1,Form("Got %d tracks on layer %s",ntTot,lr->GetName()));
+ // lr->GetMCTracks()->Print();
+ //
+ } // end loop over layers
+ //
+ // do we use vertex constraint?
+ KMCLayer *vtx = GetLayer(0);
+ if (!vtx->IsDead() && vtx->IsITS()) {
+ int ntr = vtx->GetNMCTracks();
+ for (int itr=0;itr<ntr;itr++) {
+ currTr = vtx->GetMCTrack(itr);
+ if (currTr->IsKilled()) continue;
+ KMCCluster* clv = vtx->GetMCCluster();
+ double meas[2] = {clv->GetY(),clv->GetZ()};
+ double measErr2[3] = {vtx->fPhiRes*vtx->fPhiRes,0,vtx->fZRes*vtx->fZRes};
+ double chi2v = currTr->GetPredictedChi2(meas,measErr2);
+ currTr->AddHit(vtx->GetActiveID(), chi2v, -1);
+ currTr->SetInnerLrChecked(vtx->GetActiveID());
+ if (NeedToKill(currTr)) currTr->Kill();
+ // if (vtx->IsITS()) {if (!UpdateTrack(currTr, vtx, vtx->GetMCCluster(), kFALSE)) {currTr->Kill();continue;}}
+ }
+ }
+ EliminateUnrelated();
+
+ return kTRUE;
+}
+
+//____________________________________________________________________________
+Bool_t KMCDetector::PropagateToLayer(KMCProbe* trc, KMCLayer* lr, int dir) const
+{
+ // bring the track to layer and rotat to frame normal to its surface
+ if (!trc->PropagateToR(lr->fR,fBFieldG, dir)) return kFALSE;
+ //
+ // rotate to frame with X axis normal to the surface (defined by ideal track)
+ if (!lr->IsVertex()) {
+ double pos[3];
+ trc->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;
+ if (!trc->Rotate(phi)) {
+ AliDebug(2,Form("Failed to rotate to the frame (phi:%+.3f) of layer at %.2f at XYZ: %+.3f %+.3f %+.3f",
+ phi,lr->fR,pos[0],pos[1],pos[2]));
+ if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
+ return kFALSE;
+ }
+ }
+ //
+ return kTRUE;
+}
+
+//____________________________________________________________________________
+Bool_t KMCDetector::UpdateTrack(KMCProbe* trc, KMCLayer* lr, KMCCluster* cl, Bool_t goToCluster) const
+{
+ // update track with measured cluster
+ // propagate to cluster
+ double meas[2] = {cl->GetY(),cl->GetZ()}; // ideal cluster coordinate
+ double measErr2[3] = {lr->fPhiRes*lr->fPhiRes,0,lr->fZRes*lr->fZRes};
+ //
+ if (goToCluster) if (!trc->PropagateToCluster(cl,fBFieldG)) return kFALSE; // track was not propagated to cluster frame
+ //
+ double chi2 = trc->GetPredictedChi2(meas,measErr2);
+ // if (chi2>fMaxChi2Cl) return kTRUE; // chi2 is too large
+ //
+ if (!trc->Update(meas,measErr2)) {
+ AliDebug(2,Form("layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
+ lr->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
+ if (AliLog::GetGlobalDebugLevel()>1) trc->Print("l");
+ return kFALSE;
+ }
+ trc->AddHit(lr->GetActiveID(), chi2);
+ //
+ return kTRUE;
+}
+
+//____________________________________________________________________________
+Int_t KMCDetector::GenBgClusters(KMCLayer* lr)
+{
+ // Generate fake clusters in precalculated RPhi,Z range
+ if (fNBgLimits<1) return 0; // limits were not set - no track was prolongated
+ //
+ // Fix search limits to avoid seeds which will anyway point very far from the vertex
+ double tolY = TMath::Sqrt(lr->fSig2EstD)*fMaxChi2ClSQ;
+ double tolZ = TMath::Sqrt(lr->fSig2EstZ)*fMaxChi2ClSQ;
+
+ // printf("Before: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
+ if (fBgYMin < lr->fClCorr.fY-tolY) fBgYMin = lr->fClCorr.fY-tolY;
+ if (fBgYMax > lr->fClCorr.fY+tolY) fBgYMax = lr->fClCorr.fY+tolY;
+ if (fBgZMin < lr->fClCorr.fZ-tolZ) fBgZMin = lr->fClCorr.fZ-tolZ;
+ if (fBgZMax > lr->fClCorr.fZ+tolZ) fBgZMax = lr->fClCorr.fZ+tolZ;
+ //printf("After: Y: %+6.3f : %+6.3f tolY: %6.3f || Z: %+6.3f : %+6.3f tolZ: %6.3f\n",fBgYMin,fBgYMax,tolY, fBgZMin,fBgZMax,tolZ);
+ //
+ double dy = fBgYMax - fBgYMin;
+ double dz = fBgZMax - fBgZMin;
+ double surf = dy*dz; // surface of generation
+ if (surf<0) return 0;
+ double poissProb = surf*HitDensity(lr->fR)*lr->GetLayerEff();
+ AliDebug(2,Form("Bg for Lr %s (r=%.2f) : Density %.2f on surface %.2e [%+.4f : %+.4f][%+.4f %+.4f]",
+ lr->GetName(),lr->fR,HitDensity(lr->fR),surf,fBgYMin,fBgYMax,fBgZMin,fBgZMax));
+ int nFakesGen = gRandom->Poisson( poissProb ); // preliminary number of extra clusters to test
+ KMCCluster *refCl = lr->GetMCCluster();
+ double sig2y = lr->GetPhiRes()*lr->GetPhiRes();
+ double sig2z = lr->GetZRes()*lr->GetZRes();
+ for (int ic=nFakesGen;ic--;) {
+ double y = fBgYMin+dy*gRandom->Rndm();
+ double z = fBgZMin+dz*gRandom->Rndm();
+ double dfy = y-refCl->GetY();
+ double dfz = z-refCl->GetZ();
+ double dist = (dfy*dfy)/sig2y + (dfz*dfz)/sig2z;
+ if (dist<4) continue; // avoid overlap with MC cluster
+ lr->AddBgCluster(y, z, refCl->GetX(), refCl->GetPhi());
+ }
+ AliDebug(2,Form("Added %6d noise clusters on lr %s (poisson Prob=%8.2f for surface %.2e) DY:%7.4f DZ: %7.4f",
+ lr->GetNBgClusters(),lr->GetName(),poissProb,surf,dy,dz));
+ return nFakesGen;
+ //
+}
+
+//____________________________________________________________________________
+void KMCDetector::UpdateSearchLimits(KMCProbe* probe, KMCLayer* lr)
+{
+ // define the search window for track on layer (where the bg hist will be generated)
+ static double *currYMin = fBgYMinTr.GetArray();
+ static double *currYMax = fBgYMaxTr.GetArray();
+ static double *currZMin = fBgZMinTr.GetArray();
+ static double *currZMax = fBgZMaxTr.GetArray();
+ //
+ double sizeY = probe->GetSigmaY2(), sizeZ = probe->GetSigmaZ2();
+ // /*
+ if (probe->GetNITSHits()<3) sizeY = 10*lr->fSig2EstD;
+ if (probe->GetNITSHits()<2) sizeZ = 10*lr->fSig2EstZ;
+ sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY+lr->fPhiRes*lr->fPhiRes); // max deviation in rphi to accept
+ sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ+lr->fZRes*lr->fZRes); // max deviation in dz to accept
+ // */
+ //
+ /*
+ if (probe->GetNITSHits()<3) sizeY = 1./(1./lr->fSig2EstD + 1./sizeY);
+ if (probe->GetNITSHits()<2) sizeZ = 1./(1./lr->fSig2EstZ + 1./sizeZ);
+ sizeY = fMaxChi2ClSQ*TMath::Sqrt(sizeY); // max deviation in rphi to accept
+ sizeZ = fMaxChi2ClSQ*TMath::Sqrt(sizeZ); // max deviation in dz to accept
+ */
+ //
+ // if (sizeY>2) sizeY=2;
+ // if (sizeZ>2) sizeZ=2;
+ // printf("Sizes at %s: %.5f %.5f\n",lr->GetName(), sizeY,sizeZ);
+ //
+ if (fNBgLimits>=fBgYMinTr.GetSize()) { // expand arrays, update pointers
+ fBgYMinTr.Set(2*(fNBgLimits+1));
+ fBgYMaxTr.Set(2*(fNBgLimits+1));
+ fBgZMinTr.Set(2*(fNBgLimits+1));
+ fBgZMaxTr.Set(2*(fNBgLimits+1));
+ currYMin = fBgYMinTr.GetArray();
+ currYMax = fBgYMaxTr.GetArray();
+ currZMin = fBgZMinTr.GetArray();
+ currZMax = fBgZMaxTr.GetArray();
+ }
+ if (fBgYMin > (currYMin[fNBgLimits]=probe->GetY()-sizeY) ) fBgYMin = currYMin[fNBgLimits];
+ if (fBgYMax < (currYMax[fNBgLimits]=probe->GetY()+sizeY) ) fBgYMax = currYMax[fNBgLimits];
+ if (fBgZMin > (currZMin[fNBgLimits]=probe->GetZ()-sizeZ) ) fBgZMin = currZMin[fNBgLimits];
+ if (fBgZMax < (currZMax[fNBgLimits]=probe->GetZ()+sizeZ) ) fBgZMax = currZMax[fNBgLimits];
+ if (AliLog::GetGlobalDebugLevel()>=2) {
+ probe->Print("l");
+ AliInfo(Form("Seed%3d Lr %s limits for y:%+8.4f z:%+8.4f [%+.4f : %+.4f][%+.4f %+.4f]",fNBgLimits,lr->GetName(),probe->GetY(),probe->GetZ(),currYMin[fNBgLimits],currYMax[fNBgLimits],currZMin[fNBgLimits],currZMax[fNBgLimits]));
+ AliInfo(Form("Global Limits Lr %s [%+.4f : %+.4f][%+.4f %+.4f]",lr->GetName(),fBgYMin,fBgYMax,fBgZMin,fBgZMax));
+ AliInfo(Form("MC Cluster: %+.4f : %+.4f",lr->fClMC.fY, lr->fClMC.fZ));
+ }
+ probe->SetUniqueID(fNBgLimits++);
+ //
+ if (lr->IsITS() && probe->GetNFakeITSHits()==0) {
+ if (fHMCLrResidRPhi) fHMCLrResidRPhi->Fill(probe->GetY() - lr->GetMCCluster()->GetY(), lr->GetActiveID());
+ if (fHMCLrResidZ) fHMCLrResidZ->Fill(probe->GetZ() - lr->GetMCCluster()->GetZ(),lr->GetActiveID());
+ }
+ //
+}
+
+//____________________________________________________________________________
+void KMCDetector::CheckTrackProlongations(KMCProbe *probe, KMCLayer* lr, KMCLayer* lrP)
+{
+ // explore prolongation of probe from lrP to lr with all possible clusters of lrP
+ // the probe is already brought to clusters frame
+ int nCl = lrP->GetNBgClusters();
+ double measErr2[3] = {lrP->fPhiRes*lrP->fPhiRes,0,lrP->fZRes*lrP->fZRes};
+ double meas[2] = {0,0};
+ UInt_t tmpID = probe->GetUniqueID();
+ double yMin = fBgYMinTr[tmpID];
+ double yMax = fBgYMaxTr[tmpID];
+ double zMin = fBgZMinTr[tmpID];
+ double zMax = fBgZMaxTr[tmpID];
+ //
+ probe->SetInnerLrChecked(lrP->GetActiveID());
+ for (int icl=-1;icl<nCl;icl++) {
+ KMCCluster* cl = icl<0 ? lrP->GetMCCluster() : lrP->GetBgCluster(icl); // -1 is for true MC cluster
+ if (cl->IsKilled()) {
+ if (AliLog::GetGlobalDebugLevel()>1) {printf("Skip cluster %d ",icl); cl->Print();}
+ continue;
+ }
+ double y = cl->GetY();
+ double z = cl->GetZ();
+ AliDebug(2,Form("Check seed%d against cl#%d out of %d at layer %s | y:%+8.4f z:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f]",tmpID,icl,nCl,lrP->GetName(),y,z,yMin,yMax,zMin,zMax));
+ if (AliLog::GetGlobalDebugLevel()>0) {
+ if (icl==-1 && probe->GetNFakeITSHits()==0) {
+ meas[0] = y; meas[1] = z;
+ double chi2a = probe->GetPredictedChi2(meas,measErr2);
+ if (chi2a>fMaxChi2Cl || (y<yMin || y>yMax) || (z<zMin || z>zMax)) {
+ probe->Print();
+ printf("Loosing good point (y:%+8.4f z:%+8.4f) on lr %s: chi2: %.2f | dy:%+8.4f dz:%+8.4f [%+.4f:%+.4f] [%+.4f:%+.4f] |x: %.2f %.2f | phi: %.2f %.2f\n",
+ y,z,lrP->GetName(),chi2a,y-probe->GetY(),z-probe->GetZ(),yMin,yMax,zMin,zMax, probe->GetX(), cl->GetX(), probe->GetAlpha(), cl->GetPhi());
+ }
+ }
+ }
+ if (y<yMin || y>yMax) continue; // preliminary check on Y
+ if (z<zMin || z>zMax) continue; // preliminary check on Z
+ meas[0] = y; meas[1] = z;
+ double chi2 = probe->GetPredictedChi2(meas,measErr2);
+ if (fHMCLrChi2 && probe->GetNFakeITSHits()==0 && icl==-1) fHMCLrChi2->Fill(chi2,lrP->GetActiveID());
+ AliDebug(2,Form("Seed-to-cluster chi2 = Chi2=%.2f",chi2));
+ if (chi2>fMaxChi2Cl) continue;
+ //
+ // update track copy
+ KMCProbe* newTr = lr->AddMCTrack( probe );
+ fUpdCalls++;
+ if (!newTr->Update(meas,measErr2)) {
+ AliDebug(2,Form("Layer %s: Failed to update the track by measurement {%.3f,%3f} err {%.3e %.3e %.3e}",
+ lrP->GetName(),meas[0],meas[1], measErr2[0],measErr2[1],measErr2[2]));
+ if (AliLog::GetGlobalDebugLevel()>1) newTr->Print("l");
+ newTr->Kill();
+ continue;
+ }
+ newTr->AddHit(lrP->GetActiveID(), chi2, icl);
+ if (AliLog::GetGlobalDebugLevel()>1) {
+ AliInfo("Cloned updated track is:");
+ newTr->Print();
+ }
+ if (NeedToKill(newTr)) newTr->Kill();
+ }
+ //
+}
+
+//____________________________________________________________________________
+void KMCDetector::ResetMCTracks(Int_t maxLr)
+{
+ int nl = GetNLayers();
+ if (maxLr<0 || maxLr>=nl) maxLr = nl-1;
+ for (int i=maxLr+1;i--;) GetLayer(i)->ResetMCTracks();
+}
+
+//____________________________________________________________________________
+Bool_t KMCDetector::NeedToKill(KMCProbe* probe) const
+{
+ // check if the seed at given layer (last one where update was tried)
+ // still has chances to be reconstructed
+ const Bool_t kModeKillMiss = kFALSE;
+ //
+ Bool_t kill = kFALSE;
+ while (1) {
+ int il = probe->GetInnerLayerChecked();
+ int nITS = probe->GetNITSHits();
+ int nITSMax = nITS + il; // maximum it can have
+ if (nITSMax<fMinITSHits) {
+ kill = kTRUE;
+ break;
+ } // has no chance to collect enough ITS hits
+ //
+ int ngr = fPattITS.GetSize();
+ if (ngr>0) { // check pattern
+ UInt_t patt = probe->GetHitsPatt();
+ // complete the layers not checked yet
+ for (int i=il;i--;) patt |= (0x1<<i);
+ for (int ig=ngr;ig--;)
+ if (!(((UInt_t)fPattITS[ig]) & patt)) {
+ kill = kTRUE;
+ break;
+ }
+ //
+ }
+ //
+ if (nITS>2) { // check if smallest possible norm chi2/ndf is acceptable
+ double chi2min = probe->GetChi2();
+ if (kModeKillMiss) {
+ int nMiss = fNActiveITSLayers - probe->GetInnerLayerChecked() - nITS; // layers already missed
+ chi2min = nMiss*probe->GetMissingHitPenalty();
+ }
+ chi2min /= ((nITSMax<<1)-KMCProbe::kNDOF);
+ if (chi2min>fMaxNormChi2NDF) {
+ kill = kTRUE;
+ break;
+ }
+ }
+ //
+ // loose vertex constraint
+ double dst;
+ if (nITS>=2) {
+ probe->GetZAt(0,fBFieldG,dst);
+ //printf("Zd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
+ if (TMath::Abs(dst)>10.) {
+ kill = kTRUE;
+ break;
+ }
+ }
+ if (nITS>=3) {
+ probe->GetYAt(0,fBFieldG,dst);
+ //printf("Dd (F%d): %f\n",probe->GetNFakeITSHits(),dst);
+ if (TMath::Abs(dst)>10.) {
+ kill = kTRUE;
+ break;
+ }
+ }
+ //
+ break;
+ }
+ if (kill && AliLog::GetGlobalDebugLevel()>1 && probe->GetNFakeITSHits()==0) {
+ printf("Killing good seed, last upd layer was %d\n",probe->GetInnerLayerChecked());
+ probe->Print("l");
+ }
+ return kill;
+}
+
+//_____________________________________________________________________
+void KMCDetector::EliminateUnrelated()
+{
+ // kill useless tracks
+ KMCLayer* lr = GetLayer(0);
+ int ntr = lr->GetNMCTracks();
+ int nval = 0;
+ for (int itr=0;itr<ntr;itr++) {
+ KMCProbe* probe = lr->GetMCTrack(itr);
+ if (probe->IsKilled()) continue;
+ if (probe->GetNITSHits()-probe->GetNFakeITSHits()<1) {probe->Kill(); continue;}
+ nval++;
+ }
+ lr->GetMCTracks()->Sort();
+ const int kDump = 0;
+ if (kDump>0) {
+ printf("Valid %d out of %d\n",nval, ntr);
+ ntr = ntr>kDump ? kDump:0;
+ for (int itr=0;itr<ntr;itr++) {
+ lr->GetMCTrack(itr)->Print();
+ }
+ }
+}
+
+//_____________________________________________________________________
+void KMCDetector::RequirePattern(UInt_t *patt, int groups)
+{
+ if (groups<1) {fPattITS.Set(0); return;}
+ fPattITS.Set(groups);
+ for (int i=0;i<groups;i++) fPattITS[i] = patt[i];
+}
+
+//_____________________________________________________________________
+void KMCDetector::CalcHardSearchLimits(double dzv)
+{
+ //
+ TArrayD zlims;
+ zlims.Set(fNActiveITSLayers);
+ for (int il0=0;il0<fNActiveITSLayers;il0++) {
+ KMCLayer* lr0 = GetActiveLayer(il0);
+ double angZTol = dzv/lr0->GetRadius();
+ for (int il1=0;il1<fNActiveITSLayers;il1++) {
+ if (il1==il0) continue;
+ KMCLayer* lr1 = GetActiveLayer(il1);
+ double ztol = angZTol*TMath::Abs(lr0->GetRadius() - lr1->GetRadius());
+ if (ztol>zlims[il1]) zlims[il1] = ztol;
+ }
+ }
+ //
+ for (int il=0;il<fNActiveITSLayers;il++) printf("ZTol%d: %8.4f\n",il,zlims[il]);
+}
+
+//_______________________________________________________
+double KMCDetector::PropagateBack(KMCProbe* trc)
+{
+ static KMCProbe bwd;
+ bwd = *trc;
+ bwd.ResetCovMat();
+ static double measErr2[3] = {0,0,0};
+ static double meas[2] = {0,0};
+ int icl = 0;
+ double chi2Tot = 0;
+ for (int il=1;il<=fLastActiveITSLayer;il++) {
+ KMCLayer* lr = GetLayer(il);
+ if (!PropagateToLayer(&bwd,lr,1)) return -1;
+ int aID = lr->GetActiveID();
+ if (aID>-1 && (icl=bwd.fClID[aID])>=-1) {
+ KMCCluster* clMC = icl<0 ? lr->GetMCCluster() : lr->GetBgCluster(icl);
+ if (!bwd.PropagateToCluster(clMC,fBFieldG)) return -1;
+ meas[0] = clMC->GetY(); meas[1] = clMC->GetZ();
+ measErr2[0] = lr->fPhiRes*lr->fPhiRes;
+ measErr2[2] = lr->fZRes*lr->fZRes;
+ double chi2a = bwd.GetPredictedChi2(meas,measErr2);
+ chi2Tot += chi2a;
+ printf("Chis %d (cl%+3d): t2c: %6.3f tot: %6.3f\n",aID,icl,chi2a, chi2Tot);
+ bwd.Update(meas,measErr2);
+ bwd.AddHit(aID, chi2a, icl);
+ }
+ if (!bwd.CorrectForMeanMaterial(lr,kFALSE)) return -1;
+ }
+ return chi2Tot;
+}
+
+
+//////////////////////////////////////////////////////////////////////////////////////////////
+//--------------------------------------------------------------------------------------------
+ClassImp(KMCTrackSummary)
+
+Int_t KMCTrackSummary::fgSumCounter = 0;
+
+//____________________________________________________________________________
+KMCTrackSummary::KMCTrackSummary(const char* name,const char* title, int nlr) :
+ TNamed(name,name), fNITSLayers(nlr),
+ fPattITS(0),fPattITSFakeExcl(0),fPattITSCorr(0),
+ fHMCChi2(0),fHMCSigDCARPhi(0),fHMCSigDCAZ(0),fHMCSigPt(0),fHMCNCl(0),fHMCFakePatt(0),
+ fCountAcc(0),fCountTot(0),fUpdCalls(0),
+ fRefProbe(),fAnProbe()
+{
+ // create summary structure for specific track conditions
+ //
+ SetMinMaxClITS();
+ SetMinMaxClITSFake();
+ SetMinMaxClITSCorr();
+ //
+ if (name==0 && title==0 && nlr==0) return; // default dummy contructor
+ //
+ enum {kNBinRes=1000};
+ const double kMaxChi2(10),kMaxResRPhi=0.1, kMaxResZ=0.1,kMaxResPt=0.3;
+ //
+ TString strN = name, strT = title;
+ if (strN.IsNull()) strN = Form("TrSum%d",fgSumCounter);
+ if (strT.IsNull()) strT = strN;
+ fgSumCounter++;
+ //
+ if (fNITSLayers<1) {
+ fNITSLayers=KMCProbe::GetNITSLayers();
+ if (fNITSLayers<1) {AliError("N ITS layer is not provided and not available from KMCProbe::GetNITSLayers()"); exit(1);}
+ AliInfo(Form("%s No N Layers provided, set %d",strN.Data(),fNITSLayers));
+ }
+ nlr = fNITSLayers;
+ //
+ TString nam,tit;
+ //
+ nam = Form("%s_mc_chi2",strN.Data());
+ tit = Form("%s MC #chi^{2}",strT.Data());
+ fHMCChi2 = new TH1F(nam.Data(),tit.Data(),kNBinRes,0,kMaxChi2);
+ fHMCChi2->GetXaxis()->SetTitle("#chi^{2}");
+ fHMCChi2->Sumw2();
+ //
+ nam = Form("%s_mc_DCArphi",strN.Data());
+ tit = Form("%s MC #DeltaDCA r#phi",strT.Data());
+ fHMCSigDCARPhi = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResRPhi,kMaxResRPhi);
+ fHMCSigDCARPhi->GetXaxis()->SetTitle("#Delta r#phi");
+ fHMCSigDCARPhi->Sumw2();
+ //
+ nam = Form("%s_mc_DCAz",strN.Data());
+ tit = Form("%s MC #DeltaDCA Z",strT.Data());
+ fHMCSigDCAZ = new TH1F(nam.Data(),tit.Data(),kNBinRes,-kMaxResZ,kMaxResZ);
+ fHMCSigDCAZ->GetXaxis()->SetTitle("#Delta Z");
+ fHMCSigDCAZ->Sumw2();
+ //
+ nam = Form("%s_mc_pt",strN.Data());
+ tit = Form("%s MC $Deltap_{T}/p_{T}",strT.Data());
+ fHMCSigPt = new TH1F(nam.Data(),tit.Data(),2*kNBinRes,-kMaxResPt,kMaxResPt);
+ fHMCSigPt->GetXaxis()->SetTitle("#Deltap_{T}/p_{T}");
+ fHMCSigPt->Sumw2();
+ //
+ nam = Form("%s_n_cl",strN.Data());
+ tit = Form("%s N Clusters",strT.Data());
+ fHMCNCl = new TH2F(nam.Data(),tit.Data(),nlr,1,nlr+1,nlr,1,nlr+1);
+ fHMCNCl->GetXaxis()->SetTitle("N Clusters Tot");
+ fHMCNCl->GetYaxis()->SetTitle("N Clusters Fake");
+ fHMCNCl->Sumw2();
+ //
+ nam = Form("%s_fake_cl",strN.Data());
+ tit = Form("%s Fake Clusters",strT.Data());
+ fHMCFakePatt = new TH1F(nam.Data(),tit.Data(),nlr,-0.5,nlr-0.5);
+ fHMCFakePatt->GetXaxis()->SetTitle("Fake clusters pattern");
+ fHMCFakePatt->Sumw2();
+ //
+}
+
+//____________________________________________________________________________
+KMCTrackSummary::~KMCTrackSummary()
+{
+ if (fHMCChi2) delete fHMCChi2;
+ if (fHMCSigDCARPhi) delete fHMCSigDCARPhi;
+ if (fHMCSigDCAZ) delete fHMCSigDCAZ;
+ if (fHMCSigPt) delete fHMCSigPt;
+ if (fHMCNCl) delete fHMCNCl;
+ if (fHMCFakePatt) delete fHMCFakePatt;
+}
+
+
+//____________________________________________________________________________
+void KMCTrackSummary::Add(const KMCTrackSummary* src, double scl)
+{
+ if (fHMCChi2) fHMCChi2->Add(src->fHMCChi2,scl);
+ if (fHMCSigDCARPhi) fHMCSigDCARPhi->Add(src->fHMCSigDCARPhi,scl);
+ if (fHMCSigDCAZ) fHMCSigDCAZ->Add(src->fHMCSigDCAZ,scl);
+ if (fHMCSigPt) fHMCSigPt->Add(src->fHMCSigPt,scl);
+ if (fHMCNCl) fHMCNCl->Add(src->fHMCNCl,scl);
+ if (fHMCFakePatt) fHMCFakePatt->Add(src->fHMCFakePatt,scl);
+ fCountAcc += src->fCountAcc;
+ fUpdCalls += src->fUpdCalls;
+}
+
+//____________________________________________________________________________
+void KMCTrackSummary::PrependTNamed(TNamed* obj, const char* nm, const char* tit)
+{
+ // prepend name, title by this prefix
+ TString name = obj->GetName(),title = obj->GetTitle();
+ name.Prepend(nm); title.Prepend(tit);
+ obj->SetNameTitle(name.Data(),title.Data());
+ //
+}
+
+//____________________________________________________________________________
+void KMCTrackSummary::SetNamePrefix(const char* pref)
+{
+ // prepend all names by this prefix
+ TString nmT,nmP = pref;
+ if (nmP.IsNull()) return;
+ nmP = Form("%s_",pref);
+ nmT = Form("%s ",pref);
+ PrependTNamed(this, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCChi2, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCSigDCARPhi, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCSigDCAZ, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCSigPt, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCNCl, nmP.Data(), nmT.Data());
+ PrependTNamed(fHMCFakePatt, nmP.Data(), nmT.Data());
+ //
+}
+
+//____________________________________________________________________________
+Bool_t KMCTrackSummary::CheckTrack(KMCProbe* trc)
+{
+ // check if the track satisfies to selection conditions
+ fCountTot++;
+ if (!trc) return kFALSE;
+ if (OutOfRange(trc->GetNITSHits(), fMinMaxClITS)) return kFALSE;
+ if (OutOfRange(trc->GetNFakeITSHits(), fMinMaxClITSFake)) return kFALSE;
+ if (OutOfRange(trc->GetNITSHits()-trc->GetNFakeITSHits(), fMinMaxClITSCorr)) return kFALSE;
+ //
+ // check layer patterns if requested
+ UInt_t patt = trc->GetHitsPatt();
+ UInt_t pattF = trc->GetFakesPatt();
+ UInt_t pattC = patt&(~pattF); // correct hits
+ // is there at least one hit in each of requested layer groups?
+ for (int ip=fPattITS.GetSize();ip--;) if (!CheckPattern(patt,fPattITS[ip])) return kFALSE;
+ // is there at least one fake in any of requested layer groups?
+ for (int ip=fPattITSFakeExcl.GetSize();ip--;) if (CheckPattern(pattF,fPattITSFakeExcl[ip])) return kFALSE;
+ // is there at least one correct hit in each of requested layer groups?
+ for (int ip=fPattITSCorr.GetSize();ip--;) if (!CheckPattern(pattC,fPattITSCorr[ip])) return kFALSE;
+ //
+ fCountAcc++;
+ return kTRUE;
+}
+
+//____________________________________________________________________________
+void KMCTrackSummary::AddTrack(KMCProbe* trc)
+{
+ // fill track info
+ if (!CheckTrack(trc)) return;
+ //
+ fHMCChi2->Fill(trc->GetNormChi2(kFALSE));
+ fHMCSigDCARPhi->Fill(trc->GetY());
+ fHMCSigDCAZ->Fill(trc->GetZ());
+ if (fRefProbe.Pt()>0) fHMCSigPt->Fill( trc->Pt()/fRefProbe.Pt()-1.);
+ // printf("Pts: %.3f %.3f -> %.3f\n",trc->Pt(),fRefProbe.Pt(),trc->Pt()/fRefProbe.Pt()-1.);
+ // trc->Print("l");
+ // fRefProbe.Print("l");
+ fHMCNCl->Fill(trc->GetNITSHits(),trc->GetNFakeITSHits());
+ for (int i=trc->GetNITSLayers();i--;) if (trc->IsHitFake(i)) fHMCFakePatt->Fill(i);
+ //
+}
+
+//____________________________________________________________________________
+UInt_t KMCTrackSummary::Bits(Bool_t l0, Bool_t l1, Bool_t l2, Bool_t l3, Bool_t l4, Bool_t l5, Bool_t l6, Bool_t l7,Bool_t l8, Bool_t l9,
+ Bool_t l10,Bool_t l11,Bool_t l12,Bool_t l13,Bool_t l14,Bool_t l15,Bool_t l16,Bool_t l17,Bool_t l18,Bool_t l19,
+ Bool_t l20,Bool_t l21,Bool_t l22,Bool_t l23,Bool_t l24,Bool_t l25,Bool_t l26,Bool_t l27,Bool_t l28,Bool_t l29,
+ Bool_t l30,Bool_t l31)
+{
+ // create corresponding bit pattern
+ UInt_t patt = 0;
+ if (l0 ) patt |= (0x1<<0); if (l1 ) patt |= (0x1<<1); if (l2 ) patt |= (0x1<<2);
+ if (l3 ) patt |= (0x1<<3); if (l4 ) patt |= (0x1<<4); if (l5 ) patt |= (0x1<<5);
+ if (l6 ) patt |= (0x1<<6); if (l7 ) patt |= (0x1<<7); if (l8 ) patt |= (0x1<<8);
+ if (l9 ) patt |= (0x1<<9); if (l10) patt |= (0x1<<10); if (l11) patt |= (0x1<<11);
+ if (l12) patt |= (0x1<<12); if (l13) patt |= (0x1<<13); if (l14) patt |= (0x1<<14);
+ if (l15) patt |= (0x1<<15); if (l16) patt |= (0x1<<16); if (l17) patt |= (0x1<<17);
+ if (l18) patt |= (0x1<<18); if (l19) patt |= (0x1<<19); if (l20) patt |= (0x1<<20);
+ if (l21) patt |= (0x1<<21); if (l22) patt |= (0x1<<22); if (l23) patt |= (0x1<<23);
+ if (l24) patt |= (0x1<<24); if (l25) patt |= (0x1<<25); if (l26) patt |= (0x1<<26);
+ if (l27) patt |= (0x1<<27); if (l28) patt |= (0x1<<28); if (l29) patt |= (0x1<<29);
+ if (l30) patt |= (0x1<<30); if (l31) patt |= (0x1<<31);
+ return patt;
+}
+
+//__________________________________________________________________________
+void KMCTrackSummary::Print(Option_t* ) const
+{
+ printf("%s: summary for track M=%5.3f pT: %6.3f eta: %.2f\n", GetName(),
+ fRefProbe.GetMass(),fRefProbe.Pt(), fRefProbe.Eta());
+ printf("Cuts on NCl ITS: Tot: %2d - %2d Fake: %2d - %2d Corr: %2d - %2d\n",
+ fMinMaxClITS[0],fMinMaxClITS[1],
+ fMinMaxClITSFake[0],fMinMaxClITSFake[1],
+ fMinMaxClITSCorr[0],fMinMaxClITSCorr[1]);
+ //
+ int nlr = fNITSLayers;
+ if (fPattITS.GetSize()) {
+ printf("Require at least 1 hit in groups: ");
+ printf("Hits obligatory in groups: ");
+ for (int i=fPattITS.GetSize();i--;) {
+ UInt_t pat = (UInt_t)fPattITS[i];
+ printf("[");
+ for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
+ printf("] ");
+ }
+ printf("\n");
+ }
+ //
+ if (fPattITSFakeExcl.GetSize()) {
+ printf("Fake hit forbidden in groups : ");
+ for (int i=fPattITSFakeExcl.GetSize();i--;) {
+ UInt_t pat = (UInt_t)fPattITSFakeExcl[i];
+ printf("[");
+ for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
+ printf("] ");
+ }
+ printf("\n");
+ }
+ //
+ if (fPattITSCorr.GetSize()) {
+ printf("Correct hit obligatory in groups: ");
+ for (int i=fPattITSCorr.GetSize();i--;) {
+ UInt_t pat = (UInt_t)fPattITSCorr[i];
+ printf("[");
+ for (int j=0;j<nlr;j++) if (pat&(0x1<<j)) printf("%d ",j);
+ printf("] ");
+ }
+ printf("\n");
+ }
+ //
+ printf("Entries: Tot: %.4e Acc: %.4e -> Eff: %.4f+-%.4f %.2e KMC updates\n",fCountTot,fCountAcc,GetEff(),GetEffErr(),GetUpdCalls());
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff