| 54fc064a |
1 | #include "Detector.h" |
| 2 | #include <TMath.h> |
| 3 | #include <TMatrixD.h> |
| 4 | #include <TGraph.h> |
| 5 | #include <TAxis.h> |
| 6 | #include <TFormula.h> |
| 7 | |
| 8 | /*********************************************************** |
| 9 | |
| 10 | Fast Simulation tool for Inner Tracker Systems |
| 11 | |
| 12 | original code of using the billoir technique was developed |
| 13 | for the HFT (STAR), James H. Thomas, jhthomas@lbl.gov |
| 14 | http://rnc.lbl.gov/~jhthomas |
| 15 | |
| 16 | Changes by S. Rossegger |
| 17 | Dec. 2010 - Translation into C++ class format |
| 18 | - Adding various Setters and Getters to build the geometry |
| 19 | (based on cylinders) plus handling of the layer properties |
| 20 | |
| 21 | ***********************************************************/ |
| 22 | |
| 23 | |
| 24 | #define RIDICULOUS 99999 // A ridiculously large resolution (cm) to flag a dead detector |
| 25 | |
| 26 | #define Luminosity 1.e27 // Luminosity of the beam (LHC HI == 1.e27, RHIC II == 8.e27 ) |
| 27 | #define SigmaD 6.0 // Size of the interaction diamond (cm) (LHC = 6.0 cm) |
| 28 | #define dNdEtaMinB 950 // Multiplicity per unit Eta (AuAu MinBias = 170, Central = 700) |
| 29 | #define dNdEtaCent 2300 // Multiplicity per unit Eta (LHC at 5.5 TeV not known) |
| 30 | |
| 31 | #define CrossSectionMinB 8 // minB Cross section for event under study (PbPb MinBias ~ 8 Barns) |
| 32 | #define AcceptanceOfTpcAndSi 1//0.60 //0.35 // Assumed geometric acceptance (efficiency) of the TPC and Si detectors |
| 33 | #define UPCBackgroundMultiplier 1.0 // Increase multiplicity in detector (0.0 to 1.0 * UPCRate ) (eg 1.0) |
| 34 | #define OtherBackground 0.0 // Increase multiplicity in detector (0.0 to 1.0 * minBias) (eg 0.0) |
| 35 | #define EfficiencySearchFlag 1 // Define search method. ChiSquare = 1, Simple = 0. ChiSq is better |
| 36 | |
| 37 | #define PionMass 0.139 // Mass of the Pion |
| 38 | #define KaonMass 0.498 // Mass of the Kaon |
| 39 | #define D0Mass 1.865 // Mass of the D0 |
| 40 | |
| 41 | |
| 42 | |
| 43 | class CylLayer : public TNamed { |
| 44 | public: |
| 45 | |
| 46 | CylLayer(char *name) : TNamed(name,name) {} |
| 47 | |
| 48 | Float_t GetRadius() {return radius;} |
| 49 | Float_t GetRadL() {return radL;} |
| 50 | Float_t GetPhiRes() {return phiRes;} |
| 51 | Float_t GetZRes() {return zRes;} |
| 52 | |
| 53 | // void Print() {printf(" r=%3.1lf X0=%1.6lf sigPhi=%1.4lf sigZ=%1.4lf\n",radius,radL,phiRes,zRes); } |
| 54 | Float_t radius; Float_t radL; Float_t phiRes; Float_t zRes; |
| 55 | Bool_t isDead; |
| 56 | |
| 57 | ClassDef(CylLayer,1); |
| 58 | }; |
| 59 | |
| 60 | ClassImp(Detector) |
| 61 | Detector::Detector() |
| 62 | : TNamed("test_detector","detector"), |
| 63 | numberOfLayers(0), |
| 64 | fBField(0.5), |
| 65 | fLhcUPCscale(2.5), |
| 66 | fIntegrationTime(0.02), // in ms |
| 67 | fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle |
| 68 | fParticleMass(0.140) // Standard: pion mass |
| 69 | { |
| 70 | // |
| 71 | // default constructor |
| 72 | // |
| 73 | // fLayers = new TObjArray(); |
| 74 | } |
| 75 | |
| 76 | Detector::Detector(char *name, char *title) |
| 77 | : TNamed(name,title), |
| 78 | numberOfLayers(0), |
| 79 | fBField(0.5), |
| 80 | fLhcUPCscale(2.5), |
| 81 | fIntegrationTime(0.02), // in ms |
| 82 | fAvgRapidity(0.45), // Avg rapidity, MCS calc is a function of crossing angle |
| 83 | fParticleMass(0.140) // Standard: pion mass |
| 84 | { |
| 85 | // |
| 86 | // default constructor, that set the name and title |
| 87 | // |
| 88 | // fLayers = new TObjArray(); |
| 89 | } |
| 90 | Detector::~Detector() { // |
| 91 | // virtual destructor |
| 92 | // |
| 93 | // delete fLayers; |
| 94 | } |
| 95 | |
| 96 | void Detector::AddLayer(char *name, Float_t radius, Float_t radL, Float_t phiRes, Float_t zRes) { |
| 97 | |
| 98 | CylLayer *newLayer = (CylLayer*) fLayers.FindObject(name); |
| 99 | |
| 100 | if (!newLayer) { |
| 101 | newLayer = new CylLayer(name); |
| 102 | newLayer->radius = radius; |
| 103 | newLayer->radL = radL; |
| 104 | newLayer->phiRes = phiRes; |
| 105 | newLayer->zRes = zRes; |
| 106 | |
| 107 | if (newLayer->zRes==RIDICULOUS && newLayer->zRes==RIDICULOUS) |
| 108 | newLayer->isDead = kTRUE; |
| 109 | else |
| 110 | newLayer->isDead = kFALSE; |
| 111 | |
| 112 | if (fLayers.GetEntries()==0) |
| 113 | fLayers.Add(newLayer); |
| 114 | else { |
| 115 | |
| 116 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { |
| 117 | CylLayer *l = (CylLayer*)fLayers.At(i); |
| 118 | if (radius<l->radius) { |
| 119 | fLayers.AddBefore(l,newLayer); |
| 120 | break; |
| 121 | } |
| 122 | if (radius>l->radius && (i+1)==fLayers.GetEntries() ) { |
| 123 | // even bigger then last one |
| 124 | fLayers.Add(newLayer); |
| 125 | } |
| 126 | } |
| 127 | |
| 128 | } |
| 129 | numberOfLayers += 1; |
| 130 | |
| 131 | } else { |
| 132 | printf("Layer with the name %s does already exist\n",name); |
| 133 | } |
| 134 | |
| 135 | |
| 136 | } |
| 137 | |
| 138 | void Detector::KillLayer(char *name) { |
| 139 | // |
| 140 | // Marks layer as dead. Contribution only by Material Budget |
| 141 | // |
| 142 | |
| 143 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 144 | if (!tmp) |
| 145 | printf("Layer %s not found - cannot mark as dead\n",name); |
| 146 | else { |
| 147 | tmp->phiRes = 99999; |
| 148 | tmp->zRes = 99999; |
| 149 | tmp->isDead = kTRUE; |
| 150 | } |
| 151 | } |
| 152 | |
| 153 | void Detector::SetRadius(char *name, Float_t radius) { |
| 154 | // |
| 155 | // Set layer radius [cm] |
| 156 | // |
| 157 | |
| 158 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 159 | |
| 160 | |
| 161 | if (!tmp) { |
| 162 | printf("Layer %s not found - cannot set resolution\n",name); |
| 163 | } else { |
| 164 | |
| 165 | Float_t tmpRadL = tmp->radL; |
| 166 | Float_t tmpPhiRes = tmp->phiRes; |
| 167 | Float_t tmpZRes = tmp->zRes; |
| 168 | |
| 169 | RemoveLayer(name); // so that the ordering is correct |
| 170 | AddLayer(name,radius,tmpRadL,tmpPhiRes,tmpZRes); |
| 171 | } |
| 172 | } |
| 173 | |
| 174 | Float_t Detector::GetRadius(char *name) { |
| 175 | // |
| 176 | // Return layer radius [cm] |
| 177 | // |
| 178 | |
| 179 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 180 | if (!tmp) |
| 181 | printf("Layer %s not found - cannot set resolution\n",name); |
| 182 | else |
| 183 | return tmp->radL; |
| 184 | |
| 185 | return 0; |
| 186 | } |
| 187 | |
| 188 | void Detector::SetRadiationLength(char *name, Float_t radL) { |
| 189 | // |
| 190 | // Set layer radius [cm] |
| 191 | // |
| 192 | |
| 193 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 194 | if (!tmp) |
| 195 | printf("Layer %s not found - cannot set resolution\n",name); |
| 196 | else { |
| 197 | tmp->radL = radL; |
| 198 | } |
| 199 | } |
| 200 | |
| 201 | Float_t Detector::GetRadiationLength(char *name) { |
| 202 | // |
| 203 | // Return layer radius [cm] |
| 204 | // |
| 205 | |
| 206 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 207 | if (!tmp) |
| 208 | printf("Layer %s not found - cannot set resolution\n",name); |
| 209 | else |
| 210 | return tmp->radL; |
| 211 | |
| 212 | return 0; |
| 213 | |
| 214 | } |
| 215 | |
| 216 | void Detector::SetResolution(char *name, Float_t phiRes, Float_t zRes) { |
| 217 | // |
| 218 | // Set layer resolution in [cm] |
| 219 | // |
| 220 | |
| 221 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 222 | if (!tmp) |
| 223 | printf("Layer %s not found - cannot set resolution\n",name); |
| 224 | else { |
| 225 | tmp->phiRes = phiRes; |
| 226 | tmp->zRes = zRes; |
| 227 | if (zRes==RIDICULOUS && phiRes==RIDICULOUS) |
| 228 | tmp->isDead = kTRUE; |
| 229 | else |
| 230 | tmp->isDead = kFALSE; |
| 231 | } |
| 232 | } |
| 233 | |
| 234 | Float_t Detector::GetResolution(char *name, Int_t axis) { |
| 235 | // |
| 236 | // Return layer resolution in [cm] |
| 237 | // axis = 0: resolution in rphi |
| 238 | // axis = 1: resolution in z |
| 239 | // |
| 240 | |
| 241 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 242 | if (!tmp) |
| 243 | printf("Layer %s not found - cannot set resolution\n",name); |
| 244 | else { |
| 245 | if (axis==0) return tmp->phiRes; |
| 246 | if (axis==1) return tmp->zRes; |
| 247 | printf("error: axis must be either 0 or 1 (rphi or z axis)\n"); |
| 248 | } |
| 249 | return 0; |
| 250 | } |
| 251 | |
| 252 | void Detector::RemoveLayer(char *name) { |
| 253 | // |
| 254 | // Removes a layer from the list |
| 255 | // |
| 256 | |
| 257 | CylLayer *tmp = (CylLayer*) fLayers.FindObject(name); |
| 258 | if (!tmp) |
| 259 | printf("Layer %s not found - cannot remove it\n",name); |
| 260 | else { |
| 261 | fLayers.Remove(tmp); |
| 262 | numberOfLayers -= 1; |
| 263 | } |
| 264 | } |
| 265 | |
| 266 | |
| 267 | void Detector::PrintLayout() { |
| 268 | // |
| 269 | // Prints the detector layout |
| 270 | // |
| 271 | |
| 272 | printf("Detector %s: \"%s\"\n",GetName(),GetTitle()); |
| 273 | |
| 274 | if (fLayers.GetEntries()>0) |
| 275 | printf(" Name \t\t r [cm] \t X0 \t phi & z res [um]\n"); |
| 276 | |
| 277 | CylLayer *tmp = 0; |
| 278 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { |
| 279 | tmp = (CylLayer*)fLayers.At(i); |
| 280 | |
| 281 | // don't print all the tpc layers |
| 282 | TString name(tmp->GetName()); |
| 283 | if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) continue; |
| 284 | |
| 285 | printf("%d. %s \t %03.2f \t%1.4f\t ",i, |
| 286 | tmp->GetName(), tmp->radius, tmp->radL); |
| 287 | if (tmp->phiRes==RIDICULOUS) |
| 288 | printf(" - "); |
| 289 | else |
| 290 | printf("%3.0f ",tmp->phiRes*10000); |
| 291 | if (tmp->zRes==RIDICULOUS) |
| 292 | printf(" -\n"); |
| 293 | else |
| 294 | printf("%3.0f\n",tmp->zRes*10000); |
| 295 | } |
| 296 | } |
| 297 | |
| 298 | void Detector::AddTPC(Float_t phiResMean, Float_t zResMean, Int_t skip) { |
| 299 | // |
| 300 | // Emulates the TPC |
| 301 | // |
| 302 | // skip=1: Use every padrow, skip=2: Signal in every 2nd padrow |
| 303 | |
| 304 | // % Radiation Lengths ... Average per TPC row (i.e. total/159 ) |
| 305 | Float_t radLPerRow = 0.000036; |
| 306 | |
| 307 | Float_t tpcInnerRadialPitch = 0.75 ; // cm |
| 308 | Float_t tpcMiddleRadialPitch = 1.0 ; // cm |
| 309 | Float_t tpcOuterRadialPitch = 1.5 ; // cm |
| 310 | // Float_t tpcInnerPadWidth = 0.4 ; // cm |
| 311 | // Float_t tpcMiddlePadWidth = 0.6 ; // cm |
| 312 | // Float_t tpcOuterPadWidth = 0.6 ; // cm |
| 313 | Float_t innerRows = 63 ; |
| 314 | Float_t middleRows = 64 ; |
| 315 | Float_t outerRows = 32 ; |
| 316 | Float_t tpcRows = (innerRows + middleRows + outerRows) ; |
| 317 | Float_t rowOneRadius = 85.2 ; // cm |
| 318 | Float_t row64Radius = 135.1 ; // cm |
| 319 | Float_t row128Radius = 199.2 ; // cm |
| 320 | |
| 321 | for ( Int_t k = 0 ; k < tpcRows ; k++ ) { |
| 322 | |
| 323 | Float_t rowRadius =0; |
| 324 | if (k<innerRows) |
| 325 | rowRadius = rowOneRadius + k*tpcInnerRadialPitch ; |
| 326 | else if ( k>=innerRows && k<(innerRows+middleRows) ) |
| 327 | rowRadius = row64Radius + (k-innerRows+1)*tpcMiddleRadialPitch ; |
| 328 | else if (k>=(innerRows+middleRows) && k<tpcRows ) |
| 329 | rowRadius = row128Radius + (k-innerRows-middleRows+1)*tpcOuterRadialPitch ; |
| 330 | |
| 331 | if ( k%skip == 0 ) |
| 332 | AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow,phiResMean,zResMean); |
| 333 | else |
| 334 | AddLayer(Form("tpc_%d",k),rowRadius,radLPerRow); // non "active" row |
| 335 | |
| 336 | |
| 337 | } |
| 338 | |
| 339 | // flag as dead, although resolution is ok ... makes live easier in the prints ... ;-) |
| 340 | CylLayer *tmp = 0; |
| 341 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { |
| 342 | tmp = (CylLayer*)fLayers.At(i); |
| 343 | TString name(tmp->GetName()); |
| 344 | if (name.Contains("tpc") && (!name.Contains("tpc_0")) ) tmp->isDead=kTRUE; |
| 345 | } |
| 346 | |
| 347 | |
| 348 | |
| 349 | } |
| 350 | |
| 351 | |
| 352 | Double_t Detector::ThetaMCS ( Double_t mass, Double_t RadLength, Double_t momentum ) |
| 353 | { |
| 354 | Double_t beta = momentum / TMath::Sqrt(momentum*momentum+mass*mass) ; |
| 355 | Double_t theta = 0.0 ; // Momentum and mass in GeV |
| 356 | // if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum ); |
| 357 | if ( RadLength > 0 ) theta = 0.0136 * TMath::Sqrt(RadLength) / ( beta * momentum ) * (1+0.038*TMath::Log(RadLength)) ; |
| 358 | return (theta) ; |
| 359 | } |
| 360 | |
| 361 | |
| 362 | Double_t Detector::ProbGoodHit ( Double_t Radius, Double_t SearchRadiusRPhi, Double_t SearchRadiusZ ) |
| 363 | { |
| 364 | // Based on work by Howard Wieman: http://rnc.lbl.gov/~wieman/GhostTracks.htm |
| 365 | // and http://rnc.lbl.gov/~wieman/HitFinding2D.htm |
| 366 | // This is the probability of getting a good hit using 2D Gaussian distribution function and infinite search radius |
| 367 | Double_t Sx, Sy, GoodHit ; |
| 368 | Sx = 2 * TMath::Pi() * SearchRadiusRPhi * SearchRadiusRPhi * HitDensity(Radius) ; |
| 369 | Sy = 2 * TMath::Pi() * SearchRadiusZ * SearchRadiusZ * HitDensity(Radius) ; |
| 370 | GoodHit = TMath::Sqrt(1./((1+Sx)*(1+Sy))) ; |
| 371 | return ( GoodHit ) ; |
| 372 | } |
| 373 | |
| 374 | |
| 375 | Double_t Detector::ProbGoodChiSqHit ( Double_t Radius, Double_t SearchRadiusRPhi, Double_t SearchRadiusZ ) |
| 376 | { |
| 377 | // Based on work by Victor Perevoztchikov and Howard Wieman: http://rnc.lbl.gov/~wieman/HitFinding2DXsq.htm |
| 378 | // This is the probability of getting a good hit using a Chi**2 search on a 2D Gaussian distribution function |
| 379 | Double_t Sx, GoodHit ; |
| 380 | Sx = 2 * TMath::Pi() * SearchRadiusRPhi * SearchRadiusZ * HitDensity(Radius) ; |
| 381 | GoodHit = 1./(1+Sx) ; |
| 382 | return ( GoodHit ) ; |
| 383 | } |
| 384 | |
| 385 | |
| 386 | Double_t Detector::HitDensity ( Double_t Radius ) |
| 387 | { |
| 388 | // Background (0-1) is included via 'OtherBackground' which multiplies the minBias rate by a scale factor. |
| 389 | // UPC electrons is a temporary kludge that is based on Kai Schweda's summary of Kai Hainken's MC results |
| 390 | // See K. Hencken et al. PRC 69, 054902 (2004) and PPT slides by Kai Schweda. |
| 391 | // Note that this function assumes we are working in CM and CM**2 [not meters]. |
| 392 | // Based on work by Yan Lu 12/20/2006, all radii and densities in centimeters or cm**2. |
| 393 | |
| 394 | Double_t MaxRadiusSlowDet = 10; //? // Maximum radius for slow detectors. Fast detectors |
| 395 | // and only fast detectors reside outside this radius. |
| 396 | Double_t ArealDensity = 0 ; |
| 397 | |
| 398 | if ( Radius > MaxRadiusSlowDet ) |
| 399 | { |
| 400 | ArealDensity = OneEventHitDensity(dNdEtaCent,Radius) ; // Fast detectors see central collision density (only) |
| 401 | ArealDensity += OtherBackground*OneEventHitDensity(dNdEtaMinB,Radius) ; // Increase density due to background |
| 402 | } |
| 403 | |
| 404 | if (Radius < MaxRadiusSlowDet ) |
| 405 | { // Note that IntegratedHitDensity will always be minB one event, or more, even if integration time => zero. |
| 406 | ArealDensity = OneEventHitDensity(dNdEtaCent,Radius) |
| 407 | + IntegratedHitDensity(dNdEtaMinB,Radius) |
| 408 | + UpcHitDensity(Radius) ; |
| 409 | ArealDensity += OtherBackground*IntegratedHitDensity(dNdEtaMinB,Radius) ; |
| 410 | // Increase density due to background |
| 411 | } |
| 412 | |
| 413 | return ( ArealDensity ) ; |
| 414 | } |
| 415 | |
| 416 | |
| 417 | double Detector::OneEventHitDensity( Double_t Multiplicity, Double_t Radius ) |
| 418 | { |
| 419 | // This is for one event at the vertex. No smearing. |
| 420 | double den = Multiplicity / (2.*TMath::Pi()*Radius*Radius) ; // 2 eta ? |
| 421 | // note: surface of sphere is '4*pi*r^2' |
| 422 | // surface of cylinder is '2*pi*r* h' |
| 423 | return den ; |
| 424 | } |
| 425 | |
| 426 | |
| 427 | double Detector::IntegratedHitDensity(Double_t Multiplicity, Double_t Radius) |
| 428 | { |
| 429 | // The integral of minBias events smeared over a gaussian vertex distribution. |
| 430 | // Based on work by Yan Lu 12/20/2006, all radii in centimeters. |
| 431 | |
| 432 | Double_t ZdcHz = Luminosity * 1.e-24 * CrossSectionMinB ; |
| 433 | Double_t den = ZdcHz * fIntegrationTime/1000. * Multiplicity * Dist(0., Radius) / (2.*TMath::Pi()*Radius) ; |
| 434 | |
| 435 | // Note that we do not allow the rate*time calculation to fall below one minB event at the vertex. |
| 436 | if ( den < OneEventHitDensity(Multiplicity,Radius) ) den = OneEventHitDensity(Multiplicity,Radius) ; |
| 437 | |
| 438 | return den ; |
| 439 | } |
| 440 | |
| 441 | |
| 442 | double Detector::UpcHitDensity(Double_t Radius) |
| 443 | { |
| 444 | // QED electrons ... |
| 445 | |
| 446 | Double_t UPCelectrons ; ; |
| 447 | UPCelectrons = fLhcUPCscale * (1.23 - Radius/6.5) ; // Fit to Kai Schweda summary tables at RHIC * 'scale' for LHC |
| 448 | if ( UPCelectrons < 0 ) UPCelectrons = 0.0 ; // UPC electrons fall off quickly and don't go to large R |
| 449 | UPCelectrons *= IntegratedHitDensity(dNdEtaMinB,Radius) ; // UPCs increase Mulitiplicty ~ proportional to MinBias rate |
| 450 | UPCelectrons *= UPCBackgroundMultiplier ; // Allow for an external multiplier (eg 0-1) to turn off UPC |
| 451 | |
| 452 | return UPCelectrons ; |
| 453 | } |
| 454 | |
| 455 | |
| 456 | double Detector::Dist(double z, double r) |
| 457 | { |
| 458 | // Convolute dEta/dZ distribution with assumed Gaussian of vertex z distribution |
| 459 | // Based on work by Howard Wieman http://rnc.lbl.gov/~wieman/HitDensityMeasuredLuminosity7.htm |
| 460 | // Based on work by Yan Lu 12/20/2006, all radii and Z location in centimeters. |
| 461 | Int_t Index = 1 ; // Start weight at 1 for Simpsons rule integration |
| 462 | Int_t NSTEPS = 301 ; // NSteps must be odd for Simpson's rule to work |
| 463 | double dist = 0.0 ; |
| 464 | double dz0 = ( 4*SigmaD - (-4)*SigmaD ) / (NSTEPS-1) ; //cm |
| 465 | double z0 = 0.0 ; //cm |
| 466 | for(int i=0; i<NSTEPS; i++){ |
| 467 | if ( i == NSTEPS-1 ) Index = 1 ; |
| 468 | z0 = -4*SigmaD + i*dz0 ; |
| 469 | dist += Index * (dz0/3.) * (1/sqrt(2.*TMath::Pi())/SigmaD) * exp(-z0*z0/2./SigmaD/SigmaD) * |
| 470 | (1/sqrt((z-z0)*(z-z0) + r*r)) ; |
| 471 | if ( Index != 4 ) Index = 4; else Index = 2 ; |
| 472 | } |
| 473 | return dist; |
| 474 | } |
| 475 | |
| 476 | #define PZero 0.861 // Momentum of back to back decay particles in the CM frame |
| 477 | #define EPiZero 0.872 // Energy of the pion from a D0 decay at rest |
| 478 | #define EKZero 0.993 // Energy of the Kaon from a D0 decay at rest |
| 479 | |
| 480 | Double_t Detector::D0IntegratedEfficiency( Double_t pt, Double_t corrEfficiency[][400] ) { |
| 481 | // Math from Ron Longacre. Note hardwired energy to bin conversion for PtK and PtPi. |
| 482 | |
| 483 | Double_t Const1 = pt / D0Mass ; |
| 484 | Double_t Const2 = TMath::Sqrt(pt*pt+D0Mass*D0Mass) / D0Mass ; |
| 485 | Double_t sum, PtPi, PtK ; |
| 486 | Double_t effp, effk ; |
| 487 | |
| 488 | sum = 0.0 ; |
| 489 | for ( Int_t k = 0 ; k < 360 ; k++ ) { |
| 490 | |
| 491 | Double_t theta = k * TMath::Pi() / 180. ; |
| 492 | |
| 493 | PtPi = TMath::Sqrt( |
| 494 | PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*Const2*Const2 + |
| 495 | Const1*Const1*EPiZero*EPiZero - |
| 496 | 2*PZero*TMath::Cos(theta)*Const2*Const1*EPiZero + |
| 497 | PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta) |
| 498 | ) ; |
| 499 | |
| 500 | PtK = TMath::Sqrt( |
| 501 | PZero*PZero*TMath::Cos(theta)*TMath::Cos(theta)*Const2*Const2 + |
| 502 | Const1*Const1*EKZero*EKZero + |
| 503 | 2*PZero*TMath::Cos(theta)*Const2*Const1*EKZero + |
| 504 | PZero*PZero*TMath::Sin(theta)*TMath::Sin(theta) |
| 505 | ) ; |
| 506 | |
| 507 | // JT Test Remove 100 MeV/c in pt to simulate eta!=0 decays |
| 508 | Int_t pionindex = (int)((PtPi-0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; |
| 509 | Int_t kaonindex = (int)((PtK -0.1)*100.0 - 65.0*TMath::Abs(fBField)) ; |
| 510 | |
| 511 | if ( pionindex >= 400 ) pionindex = 399 ; |
| 512 | if ( pionindex >= 0 ) effp = corrEfficiency[0][pionindex] ; |
| 513 | if ( pionindex < 0 ) effp = (corrEfficiency[0][1]-corrEfficiency[0][0])*pionindex + corrEfficiency[0][0] ; // Extrapolate if reqd |
| 514 | if ( effp < 0 ) effp = 0 ; |
| 515 | |
| 516 | if ( kaonindex >= 400 ) kaonindex = 399 ; |
| 517 | if ( kaonindex >= 0 ) effk = corrEfficiency[1][kaonindex] ; |
| 518 | if ( kaonindex < 0 ) effk = (corrEfficiency[1][1]-corrEfficiency[1][0])*kaonindex + corrEfficiency[1][0] ; // Extrapolate if reqd |
| 519 | if ( effk < 0 ) effk = 0 ; |
| 520 | |
| 521 | // Note that we assume that the Kaon Decay efficiency has already been inlcuded in the kaon efficiency used here. |
| 522 | |
| 523 | sum += effp * effk ; |
| 524 | |
| 525 | } |
| 526 | |
| 527 | Double_t mean =sum/360; |
| 528 | |
| 529 | |
| 530 | |
| 531 | return mean ; |
| 532 | |
| 533 | } |
| 534 | |
| 535 | |
| 536 | |
| 537 | |
| 538 | void Detector::SolveViaBilloir(Int_t flagD0,Int_t print) { |
| 539 | |
| 540 | // Calculate track parameters using Billoirs method of matrices |
| 541 | |
| 542 | |
| 543 | Double_t pt, pz, lambda, momentum, rho, DeltaPoverP ; |
| 544 | Double_t layerThickness, MCS, mmm, Charge ; |
| 545 | Double_t F1, F2, F3, G1, G2, G3 ; |
| 546 | Double_t dx, sigma2, sigma2Z ; |
| 547 | Double_t mass[3] ; |
| 548 | Int_t PrintOnce = 1 ; |
| 549 | |
| 550 | mass[0] = PionMass ; mass[1] = KaonMass ; // Loop twice for the D0; first pi then k |
| 551 | |
| 552 | mass[2] = fParticleMass; // third loop |
| 553 | |
| 554 | Int_t mStart =0; |
| 555 | if (!flagD0) mStart = 2; // pion and kaon is skipped -> fast mode |
| 556 | |
| 557 | for ( Int_t massloop = mStart ; massloop < 3 ; massloop++ ) { |
| 558 | |
| 559 | // PseudoRapidity OK, used as an angle |
| 560 | lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)) ; |
| 561 | |
| 562 | // Track moves along the x axis and deviations are in the y and z directions. |
| 563 | Double_t y, z, a, b, d, e ; |
| 564 | // a = py/px = tan phi, b = pz/px = tan lambda / cos phi, d = 0.3*Ze/Abs(p) |
| 565 | |
| 566 | for ( Int_t i = 0 ; i < 400 ; i++ ) { // pt loop |
| 567 | |
| 568 | CylLayer *last = (CylLayer*) fLayers.At((fLayers.GetEntries()-1)); |
| 569 | |
| 570 | // Convert to Meters, Tesla, and GeV |
| 571 | xpoint[i] = ( 0.3*(last->radius/100)*TMath::Abs(fBField) ) |
| 572 | - 0.08 + TMath::Power(10,2.3*i/400) / 10.0 ; // Starting values based on radius in TPC ... log10 steps to ~20 GeV |
| 573 | |
| 574 | |
| 575 | |
| 576 | // Assume track started at (0,0,0) and shoots out on the X axis, and B field is on the Z axis |
| 577 | // These are the EndPoint values for y, z, a, b, and d |
| 578 | pt = xpoint[i] ; // GeV/c |
| 579 | rho = pt / TMath::Abs( 0.3 * fBField ); // Radius of curvature of the track (meters) |
| 580 | momentum = pt / TMath::Cos(lambda) ; // Total momentum |
| 581 | Charge = -1 ; // Assume an electron |
| 582 | pz = pt * TMath::Tan(lambda) ; // GeV/c |
| 583 | d = 0.3 * Charge / momentum ; // Its an electrons so q = -1, which makes d a negative number |
| 584 | |
| 585 | // printf("%d - pt %lf r%lf | %lf %lf\n",massloop,xpoint[i],(last->radius)/100,momentum, d); |
| 586 | |
| 587 | // Matrices for Billoir |
| 588 | TMatrixD BigA(2,2); TMatrixD LittleA(2,2); |
| 589 | TMatrixD Istar(5,5); TMatrixD D(5,5); TMatrixD DT(5,5); TMatrixD A(5,5); TMatrixD M(5,5); TMatrixD Work(5,5) ; |
| 590 | |
| 591 | // Set Detector-Efficiency Storage area to unity |
| 592 | efficiency[massloop][i] = 1.0 ; |
| 593 | // Back-propagate the covariance matrix along the track. I is the inverse of the covariance matrix. |
| 594 | // Start with small variations to a straight line 'information' matrix. Tuning this matrix is |
| 595 | // equivalent to 'starting the recursion' in Billoirs paper. It is a tricky business. You need test |
| 596 | // data and many layers for the matrices to stabilize. Do not believe results for a small number of layers. |
| 597 | // In our case, always include the TPC and this matrix will be well conditioned before it gets to the Si. |
| 598 | |
| 599 | Istar.UnitMatrix(); // start with unity |
| 600 | |
| 601 | // 'A' is the "angle" matrix for MCS at each step. |
| 602 | // 'D' is the "distance" matrix at each step. It propagates the particle backwards along the track. |
| 603 | // 'M' is the "measurement" matrix. It represents the measurement resolution at each step. |
| 604 | |
| 605 | CylLayer *layer = 0; |
| 606 | CylLayer *nextlayer = 0; |
| 607 | |
| 608 | for (Int_t j=(fLayers.GetEntries()-1); j>0; j--) { // Layer loop |
| 609 | |
| 610 | layer = (CylLayer*)fLayers.At(j); |
| 611 | nextlayer = (CylLayer*)fLayers.At(j-1); |
| 612 | |
| 613 | // layer->Print(); |
| 614 | // Convert to Meters, Tesla, and GeV |
| 615 | Float_t radius = layer->radius/100; |
| 616 | // Float_t phiResolution = layer->GetPhiRes()/100; |
| 617 | // Float_t zResolution = layer->GetZRes()/100; |
| 618 | Float_t radLength = layer->radL; |
| 619 | |
| 620 | Float_t nextRadius = nextlayer->radius /100; |
| 621 | Float_t nextPhiResolution = nextlayer->phiRes /100; |
| 622 | Float_t nextZResolution = nextlayer->zRes /100; |
| 623 | |
| 624 | // printf("%d: %lf %lf %lf %lf\n",j,radius,nextRadius, radLength,nextPhiResolution); |
| 625 | // printf("%d %d %d : %lf %lf\n",j,fLayers.GetEntries()-1,0, rho, radius); |
| 626 | |
| 627 | if ( radius >= rho ) // (meters) protect agains too low pt |
| 628 | { printf("pt lower bound is too low for this algorithm\n"); return ; } |
| 629 | |
| 630 | y = rho - TMath::Sqrt(rho*rho-radius*radius) ; // These are 'ideal' locations and |
| 631 | a = radius / ( rho - y ) ; // not propagated values which should |
| 632 | b = rho * TMath::Tan(lambda) / ( rho - y ) ; // be done if we had data. But we don't. |
| 633 | z = rho * TMath::Tan(lambda) * TMath::ATan(a) ; |
| 634 | e = TMath::Sqrt( 1 + a*a + b*b ) ; |
| 635 | |
| 636 | layerThickness = radLength / TMath::Sin(TMath::Pi()/2 - lambda) ; |
| 637 | |
| 638 | MCS = ThetaMCS(mass[massloop], layerThickness, momentum) ; |
| 639 | mmm = ( 1 + a*a + b*b ) ; |
| 640 | BigA(0,0) = mmm*(1+a*a) ; BigA(0,1) = mmm*a*b ; LittleA(0,0) = Istar(2,2) ; LittleA(0,1) = Istar(2,3) ; |
| 641 | BigA(1,0) = mmm*a*b ; BigA(1,1) = mmm*(1+b*b) ; LittleA(1,0) = Istar(3,2) ; LittleA(1,1) = Istar(3,3) ; |
| 642 | LittleA = BigA.Invert() + MCS*MCS*LittleA ; |
| 643 | LittleA = LittleA.Invert() ; |
| 644 | A(0,0) = 0.0 ; A(0,1) = 0.0 ; A(0,2) = 0.0 ; A(0,3) = 0.0 ; A(0,4) = 0.0 ; |
| 645 | A(1,0) = 0.0 ; A(1,1) = 0.0 ; A(1,2) = 0.0 ; A(1,3) = 0.0 ; A(1,4) = 0.0 ; |
| 646 | A(2,0) = 0.0 ; A(2,1) = 0.0 ; A(2,2) = LittleA(0,0) ; A(2,3) = LittleA(0,1) ; A(2,4) = 0.0 ; |
| 647 | A(3,0) = 0.0 ; A(3,1) = 0.0 ; A(3,2) = LittleA(1,0) ; A(3,3) = LittleA(1,1) ; A(3,4) = 0.0 ; |
| 648 | A(4,0) = 0.0 ; A(4,1) = 0.0 ; A(4,2) = 0.0 ; A(4,3) = 0.0 ; A(4,4) = 0.0 ; |
| 649 | Istar = Istar - MCS*MCS*Istar*A*Istar ; |
| 650 | dx = ( radius - nextRadius ) ; // Billoir works with absolute magnitude of distance |
| 651 | F3 = e * ( -1 * ( 1 + a*a ) * fBField ) ; // Assume fBField is on Z axis, only. |
| 652 | F1 = d * ( a*F3/(e*e) - 2*e*a*fBField ) ; |
| 653 | F2 = d * ( b*F3/(e*e) ) ; |
| 654 | G3 = e * ( -1 * a * b * fBField ) ; |
| 655 | G1 = d * ( a*G3/(e*e) + e*b*fBField ) ; |
| 656 | G2 = d * ( b*G3/(e*e) - e*a*fBField ) ; |
| 657 | D(0,0) = 1.0 ; D(0,1) = 0.0 ; D(0,2) = -1*dx+F1*dx*dx/2 ; D(0,3) = F2*dx*dx/2 ; D(0,4) = F3*dx*dx/2 ; |
| 658 | D(1,0) = 0.0 ; D(1,1) = 1.0 ; D(1,2) = G1*dx*dx/2 ; D(1,3) = -1*dx+G2*dx*dx/2 ; D(1,4) = G3*dx*dx/2 ; |
| 659 | D(2,0) = 0.0 ; D(2,1) = 0.0 ; D(2,2) = 1-F1*dx ; D(2,3) = -1*F2*dx ; D(2,4) = -1*F3*dx ; |
| 660 | D(3,0) = 0.0 ; D(3,1) = 0.0 ; D(3,2) = -1*G1*dx ; D(3,3) = 1-G2*dx ; D(3,4) = -1*G3*dx ; |
| 661 | D(4,0) = 0.0 ; D(4,1) = 0.0 ; D(4,2) = 0.0 ; D(4,3) = 0.0 ; D(4,4) = 1.0 ; |
| 662 | DT.Transpose(D) ; |
| 663 | Istar = DT.Invert() * Istar * D.Invert() ; |
| 664 | // Prepare to save Detector efficiencies |
| 665 | Work = Istar ; // Working copy of matrix |
| 666 | Work.Invert() ; // Invert the Matrix to recover the convariance matrix |
| 667 | DetPointRes [j-1][i] = TMath::Sqrt( Work(0,0) ) ; // result in meters |
| 668 | DetPointZRes[j-1][i] = TMath::Sqrt( Work(1,1) ) ; // result in meters |
| 669 | // End save |
| 670 | sigma2 = ( nextPhiResolution * nextPhiResolution ) ; |
| 671 | sigma2Z = ( nextZResolution * nextZResolution ) ; |
| 672 | M(0,0) = 1/sigma2 ; M(0,1) = 0.0 ; M(0,2) = 0.0 ; M(0,3) = 0.0 ; M(0,4) = 0.0 ; |
| 673 | M(1,0) = 0.0 ; M(1,1) = 1/sigma2Z ; M(1,2) = 0.0 ; M(1,3) = 0.0 ; M(1,4) = 0.0 ; |
| 674 | M(2,0) = 0.0 ; M(2,1) = 0.0 ; M(2,2) = 0.0 ; M(2,3) = 0.0 ; M(2,4) = 0.0 ; |
| 675 | M(3,0) = 0.0 ; M(3,1) = 0.0 ; M(3,2) = 0.0 ; M(3,3) = 0.0 ; M(3,4) = 0.0 ; |
| 676 | M(4,0) = 0.0 ; M(4,1) = 0.0 ; M(4,2) = 0.0 ; M(4,3) = 0.0 ; M(4,4) = 0.0 ; |
| 677 | Istar = Istar + M ; |
| 678 | |
| 679 | } |
| 680 | |
| 681 | // Invert the Matrix to recover the convariance matrix |
| 682 | Istar.Invert() ; |
| 683 | // Convert the Convariance matrix parameters into physical quantities |
| 684 | // The results are propogated to the previous point but *do not* include the measurement at that point. |
| 685 | DeltaPoverP = TMath::Sqrt( Istar(4,4) ) * momentum / 0.3 ; // Absolute magnitude so ignore Charge |
| 686 | ypoint[i] = 100.* TMath::Abs( DeltaPoverP ) ; // results in percent |
| 687 | yprime[i] = TMath::Sqrt( Istar(0,0) ) * 1.e6 ; // result in microns |
| 688 | yprimeZ[i] = TMath::Sqrt( Istar(1,1) ) * 1.e6 ; // result in microns |
| 689 | equivalent[i] = TMath::Sqrt(yprime[i]*yprimeZ[i]) ; // Equivalent circular radius |
| 690 | |
| 691 | |
| 692 | if (print == 1 && xpoint[i] > 0.750 && massloop == 2 && PrintOnce == 1) { |
| 693 | printf("Mass of tracked particle: %f\n",fParticleMass) ; |
| 694 | printf("Radius Thickness PointResOn PointResOnZ DetRes DetResZ Density Efficiency\n") ; |
| 695 | // PrintOnce =0; |
| 696 | } |
| 697 | |
| 698 | // print and efficiency |
| 699 | for (Int_t j=(fLayers.GetEntries()-1); j>0; j--) { // Layer loop |
| 700 | |
| 701 | layer = (CylLayer*)fLayers.At(j-1); |
| 702 | |
| 703 | // Convert to Meters, Tesla, and GeV |
| 704 | Float_t radius = layer->radius /100; |
| 705 | Float_t phiRes = layer->phiRes /100; |
| 706 | Float_t zRes = layer->zRes /100; |
| 707 | Float_t radLength = layer->radL; |
| 708 | Bool_t isDead = layer->isDead; |
| 709 | |
| 710 | if ( !isDead && radLength >0 ) { |
| 711 | Double_t rphiError = TMath::Sqrt( DetPointRes[j-1][i] * DetPointRes [j-1][i] + |
| 712 | phiRes * phiRes ) * 100. ; // work in cm |
| 713 | Double_t zError = TMath::Sqrt( DetPointZRes[j-1][i] * DetPointZRes[j-1][i] + |
| 714 | zRes * zRes ) * 100. ; // work in cm |
| 715 | |
| 716 | if (print == 1 && xpoint[i] > 0.750 && massloop == 2 && PrintOnce == 1) { |
| 717 | printf("%5.1f %9.4f %10.0f %11.0f %7.0f %8.0f %8.2f %10.3f\n", |
| 718 | radius*100, radLength, |
| 719 | DetPointRes[j-1][i]*1.e6, DetPointZRes[j-1][i]*1.e6, |
| 720 | phiRes*1.e6, zRes*1.e6, |
| 721 | HitDensity(radius*100), |
| 722 | ProbGoodChiSqHit( radius*100, rphiError , zError ) ) ; |
| 723 | } |
| 724 | |
| 725 | if ( EfficiencySearchFlag == 0 ) |
| 726 | efficiency[massloop][i] *= ProbGoodHit( radius*100, rphiError , zError ) ; |
| 727 | else if ( EfficiencySearchFlag == 1 ) |
| 728 | efficiency[massloop][i] *= ProbGoodChiSqHit( radius*100, rphiError , zError ) ; |
| 729 | } |
| 730 | } |
| 731 | if (print == 1 && xpoint[i] > 0.750 && massloop == 2 && PrintOnce == 1) { |
| 732 | PrintOnce = 0 ; |
| 733 | printf("\n") ; |
| 734 | } |
| 735 | |
| 736 | } // mass loop |
| 737 | } // pt loop |
| 738 | |
| 739 | } |
| 740 | |
| 741 | |
| 742 | TGraph * Detector::GetGraphMomentumResolution(Int_t color, Int_t linewidth) { |
| 743 | |
| 744 | TGraph *graph = new TGraph(400, xpoint, ypoint); |
| 745 | graph->SetMaximum(10) ; |
| 746 | graph->SetMinimum(0) ; |
| 747 | graph->SetTitle("Momentum Resolution .vs. Pt" ) ; |
| 748 | graph->GetXaxis()->SetRangeUser(0.,5.0) ; |
| 749 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; |
| 750 | graph->GetXaxis()->CenterTitle(); |
| 751 | graph->GetYaxis()->SetTitle("Momentum Resolution [%]") ; |
| 752 | graph->GetYaxis()->CenterTitle(); |
| 753 | |
| 754 | graph->SetLineColor(color); |
| 755 | graph->SetMarkerColor(color); |
| 756 | graph->SetLineWidth(linewidth); |
| 757 | |
| 758 | return graph; |
| 759 | |
| 760 | } |
| 761 | |
| 762 | TGraph * Detector::GetGraphPointingResolution(Int_t axis, Int_t color, Int_t linewidth) { |
| 763 | |
| 764 | // Returns the pointing resolution |
| 765 | // axis = 0 ... rphi pointing resolution |
| 766 | // axis = 1 ... z pointing resolution |
| 767 | // |
| 768 | |
| 769 | TGraph * graph = 0; |
| 770 | |
| 771 | if (axis==0) { |
| 772 | graph = new TGraph ( 400, xpoint, yprime ) ; |
| 773 | graph->SetTitle("R-#phi Pointing Resolution .vs. Pt" ) ; |
| 774 | graph->GetYaxis()->SetTitle("R-#phi Pointing Resolution [#mum]") ; |
| 775 | } else { |
| 776 | graph = new TGraph ( 400, xpoint, yprimeZ ) ; |
| 777 | graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ; |
| 778 | graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum]") ; |
| 779 | } |
| 780 | |
| 781 | graph->SetMinimum(1) ; |
| 782 | graph->SetMaximum(300.1) ; |
| 783 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; |
| 784 | graph->GetXaxis()->CenterTitle(); |
| 785 | graph->GetXaxis()->SetNoExponent(1) ; |
| 786 | graph->GetXaxis()->SetMoreLogLabels(1) ; |
| 787 | graph->GetYaxis()->CenterTitle(); |
| 788 | |
| 789 | graph->SetLineWidth(linewidth); |
| 790 | graph->SetLineColor(color); |
| 791 | graph->SetMarkerColor(color); |
| 792 | |
| 793 | return graph; |
| 794 | |
| 795 | } |
| 796 | |
| 797 | |
| 798 | TGraph * Detector::GetGraphPointingResolutionTeleEqu(Int_t axis,Int_t color, Int_t linewidth) { |
| 799 | // |
| 800 | // returns the Pointing resolution (accoring to Telescope equation) |
| 801 | // axis =0 ... in rphi |
| 802 | // axis =1 ... in z |
| 803 | // |
| 804 | |
| 805 | Double_t resolution[400]; |
| 806 | |
| 807 | Double_t layerResolution[2]; |
| 808 | Double_t layerRadius[2]; |
| 809 | Double_t layerThickness[2]; |
| 810 | |
| 811 | Int_t count =0; // search two first active layers |
| 812 | printf("Telescope equation for layers:\n"); |
| 813 | for (Int_t i = 0; i<fLayers.GetEntries(); i++) { |
| 814 | CylLayer *l = (CylLayer*)fLayers.At(i); |
| 815 | if (!l->isDead && l->radius>0) { |
| 816 | layerRadius[count] = l->radius; |
| 817 | layerThickness[count] = l->radL; |
| 818 | if (axis==0) { |
| 819 | layerResolution[count] = l->phiRes; |
| 820 | } else { |
| 821 | layerResolution[count] = l->zRes; |
| 822 | } |
| 823 | printf(" %s \n",l->GetName()); |
| 824 | count++; |
| 825 | } |
| 826 | if (count>=2) break; |
| 827 | } |
| 828 | |
| 829 | |
| 830 | Double_t pt, momentum, thickness,MCS ; |
| 831 | Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); |
| 832 | |
| 833 | for ( Int_t i = 0 ; i < 400 ; i++ ) { |
| 834 | // Reference data as if first two layers were acting all alone |
| 835 | pt = xpoint[i] ; |
| 836 | momentum = pt / TMath::Cos(lambda) ; // Total momentum |
| 837 | resolution[i] = layerResolution[0]*layerResolution[0]*layerRadius[1]*layerRadius[1] |
| 838 | + layerResolution[1]*layerResolution[1]*layerRadius[0]*layerRadius[0] ; |
| 839 | resolution[i] /= ( layerRadius[1] - layerRadius[0] ) * ( layerRadius[1] - layerRadius[0] ) ; |
| 840 | thickness = layerThickness[0] / TMath::Sin(TMath::Pi()/2 - lambda) ; |
| 841 | MCS = ThetaMCS(fParticleMass, thickness, momentum) ; |
| 842 | resolution[i] += layerRadius[0]*layerRadius[0]*MCS*MCS ; |
| 843 | resolution[i] = TMath::Sqrt(resolution[i]) * 10000.0 ; // result in microns |
| 844 | } |
| 845 | |
| 846 | |
| 847 | |
| 848 | TGraph* graph = new TGraph ( 400, xpoint, resolution ) ; |
| 849 | |
| 850 | if (axis==0) { |
| 851 | graph->SetTitle("RPhi Pointing Resolution .vs. Pt" ) ; |
| 852 | graph->GetYaxis()->SetTitle("RPhi Pointing Resolution [#mum] ") ; |
| 853 | } else { |
| 854 | graph->SetTitle("Z Pointing Resolution .vs. Pt" ) ; |
| 855 | graph->GetYaxis()->SetTitle("Z Pointing Resolution [#mum] ") ; |
| 856 | } |
| 857 | graph->SetMinimum(1) ; |
| 858 | graph->SetMaximum(300.1) ; |
| 859 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; |
| 860 | graph->GetXaxis()->CenterTitle(); |
| 861 | graph->GetXaxis()->SetNoExponent(1) ; |
| 862 | graph->GetXaxis()->SetMoreLogLabels(1) ; |
| 863 | graph->GetYaxis()->CenterTitle(); |
| 864 | |
| 865 | graph->SetLineColor(color); |
| 866 | graph->SetMarkerColor(color); |
| 867 | graph->SetLineStyle(kDashed); |
| 868 | graph->SetLineWidth(linewidth); |
| 869 | |
| 870 | return graph; |
| 871 | |
| 872 | } |
| 873 | |
| 874 | TGraph * Detector::GetGraphRecoEfficiency(Int_t particle,Int_t color, Int_t linewidth) { |
| 875 | // |
| 876 | // particle = 0 ... choosen particle (setted particleMass) |
| 877 | // particle = 1 ... Pion |
| 878 | // particle = 2 ... Kaon |
| 879 | // particle = 3 ... D0 |
| 880 | // |
| 881 | Double_t lambda = TMath::Pi()/2.0 - 2.0*TMath::ATan(TMath::Exp(-1*fAvgRapidity)); |
| 882 | |
| 883 | Double_t particleEfficiency[400]; // with chosen particle mass |
| 884 | Double_t kaonEfficiency[400], pionEfficiency[400], D0efficiency[400]; |
| 885 | Double_t partEfficiency[2][400]; |
| 886 | |
| 887 | if (particle != 0) { |
| 888 | // resulting Pion and Kaon efficiency scaled with overall efficiency |
| 889 | Double_t doNotDecayFactor; |
| 890 | for ( Int_t massloop = 0 ; massloop < 2 ; massloop++) { //0-pion, 1-kaon |
| 891 | |
| 892 | for ( Int_t j = 0 ; j < 400 ; j++ ) { |
| 893 | // JT Test Let the kaon decay. If it decays inside the TPC ... then it is gone; for all decays < 130 cm. |
| 894 | Double_t momentum = xpoint[j] / TMath::Cos(lambda) ; // Total momentum at average rapidity |
| 895 | if ( massloop == 1 ) { // KAON |
| 896 | doNotDecayFactor = TMath::Exp(-130/(371*momentum/KaonMass)) ; // Decay length for kaon is 371 cm. |
| 897 | kaonEfficiency[j] = efficiency[1][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ; |
| 898 | } else { // PION |
| 899 | doNotDecayFactor = 1.0 ; |
| 900 | pionEfficiency[j] = efficiency[0][j] * AcceptanceOfTpcAndSi*doNotDecayFactor ; |
| 901 | } |
| 902 | partEfficiency[0][j] = pionEfficiency[j]; |
| 903 | partEfficiency[1][j] = kaonEfficiency[j]; |
| 904 | } |
| 905 | } |
| 906 | |
| 907 | // resulting estimate of the D0 efficiency |
| 908 | for ( Int_t j = 0 ; j < 400 ; j++ ) { |
| 909 | D0efficiency[j] = D0IntegratedEfficiency(xpoint[j],partEfficiency); |
| 910 | } |
| 911 | } else { |
| 912 | for ( Int_t j = 0 ; j < 400 ; j++ ) { |
| 913 | particleEfficiency[j] = efficiency[2][j]* AcceptanceOfTpcAndSi; |
| 914 | // NOTE: Decay factor (see kaon) should be included to be realiable |
| 915 | } |
| 916 | } |
| 917 | |
| 918 | TGraph * graph = 0; |
| 919 | if (particle==0) { |
| 920 | graph = new TGraph ( 400, xpoint, particleEfficiency ) ; // choosen mass |
| 921 | graph->SetLineWidth(1); |
| 922 | } else if (particle==1) { |
| 923 | graph = new TGraph ( 400, xpoint, pionEfficiency ) ; |
| 924 | graph->SetLineWidth(1); |
| 925 | } else if (particle ==2) { |
| 926 | graph = new TGraph ( 400, xpoint, kaonEfficiency ) ; |
| 927 | graph->SetLineWidth(1); |
| 928 | } else if (particle ==3) { |
| 929 | graph = new TGraph ( 400, xpoint, D0efficiency ) ; |
| 930 | graph->SetLineStyle(kDashed); |
| 931 | } else |
| 932 | return 0; |
| 933 | |
| 934 | graph->GetXaxis()->SetTitle("Transverse Momentum (GeV/c)") ; |
| 935 | graph->GetXaxis()->CenterTitle(); |
| 936 | graph->GetXaxis()->SetNoExponent(1) ; |
| 937 | graph->GetXaxis()->SetMoreLogLabels(1) ; |
| 938 | graph->GetYaxis()->SetTitle("Efficiency (arbitrary units)") ; |
| 939 | graph->GetYaxis()->CenterTitle(); |
| 940 | |
| 941 | graph->SetMinimum(0.01) ; |
| 942 | graph->SetMaximum(1.0) ; |
| 943 | |
| 944 | graph->SetLineColor(color); |
| 945 | graph->SetMarkerColor(color); |
| 946 | graph->SetLineWidth(linewidth); |
| 947 | |
| 948 | return graph; |
| 949 | } |
| 950 | |
| 951 | |
| 952 | TGraph* Detector::GetGraphImpactParam(Int_t mode, Int_t axis, Int_t color, Int_t linewidth) { |
| 953 | // |
| 954 | // returns the Impact Parameter d0 (convolution of pointing resolution and vtx resolution) |
| 955 | // mode 0: impact parameter (convolution of pointing and vertex resolution) |
| 956 | // mode 1: pointing resolution |
| 957 | // mode 2: vtx resolution |
| 958 | |
| 959 | |
| 960 | TGraph *graph = new TGraph(); |
| 961 | |
| 962 | // TFormula vtxResRPhi("vtxRes","50-2*x"); // 50 microns at pt=0, 15 microns at pt =20 ? |
| 963 | TFormula vtxResRPhi("vtxRes","35/(x+1)+10"); // |
| 964 | TFormula vtxResZ("vtxResZ","600/(x+4)+10"); // |
| 965 | |
| 966 | TGraph *trackRes = GetGraphPointingResolution(axis,1); |
| 967 | Double_t *pt = trackRes->GetX(); |
| 968 | Double_t *trRes = trackRes->GetY(); |
| 969 | for (Int_t ip =0; ip<trackRes->GetN(); ip++) { |
| 970 | Double_t vtxRes = 0; |
| 971 | if (axis==0) |
| 972 | vtxRes = vtxResRPhi.Eval(pt[ip]); |
| 973 | else |
| 974 | vtxRes = vtxResZ.Eval(pt[ip]); |
| 975 | |
| 976 | if (mode==0) |
| 977 | graph->SetPoint(ip,pt[ip],TMath::Sqrt(vtxRes*vtxRes+trRes[ip]*trRes[ip])); |
| 978 | else if (mode ==1) |
| 979 | graph->SetPoint(ip,pt[ip],trRes[ip]); |
| 980 | else |
| 981 | graph->SetPoint(ip,pt[ip],vtxRes); |
| 982 | } |
| 983 | |
| 984 | graph->SetTitle("d_{0} r#phi resolution .vs. Pt" ) ; |
| 985 | graph->GetYaxis()->SetTitle("d_{0} r#phi resolution [#mum]") ; |
| 986 | |
| 987 | graph->SetMinimum(1) ; |
| 988 | graph->SetMaximum(300.1) ; |
| 989 | graph->GetXaxis()->SetTitle("Transverse Momentum [GeV/c]") ; |
| 990 | graph->GetXaxis()->CenterTitle(); |
| 991 | graph->GetXaxis()->SetNoExponent(1) ; |
| 992 | graph->GetXaxis()->SetMoreLogLabels(1) ; |
| 993 | graph->GetYaxis()->CenterTitle(); |
| 994 | |
| 995 | graph->SetLineColor(color); |
| 996 | graph->SetMarkerColor(color); |
| 997 | graph->SetLineWidth(linewidth); |
| 998 | |
| 999 | return graph; |
| 1000 | |
| 1001 | } |
| 1002 | |
| 1003 | TGraph* Detector::GetGraph(Int_t number, Int_t color, Int_t linewidth) { |
| 1004 | // |
| 1005 | // returns graph according to the number |
| 1006 | // |
| 1007 | switch(number) { |
| 1008 | case 1: |
| 1009 | return GetGraphPointingResolution(0,color, linewidth); // dr |
| 1010 | case 2: |
| 1011 | return GetGraphPointingResolution(1,color, linewidth); // dz |
| 1012 | case 3: |
| 1013 | return GetGraphPointingResolutionTeleEqu(0,color, linewidth); // dr - tele |
| 1014 | case 4: |
| 1015 | return GetGraphPointingResolutionTeleEqu(1,color, linewidth); // dz - tele |
| 1016 | case 5: |
| 1017 | return GetGraphMomentumResolution(color, linewidth); // pt resolution |
| 1018 | case 11: |
| 1019 | return GetGraphRecoEfficiency(1, color, linewidth); // eff. pion |
| 1020 | case 12: |
| 1021 | return GetGraphRecoEfficiency(2, color, linewidth); // eff. kaon |
| 1022 | case 13: |
| 1023 | return GetGraphRecoEfficiency(3, color, linewidth); // eff. D0 |
| 1024 | default: |
| 1025 | printf(" Error: chosen graph number not valid\n"); |
| 1026 | } |
| 1027 | return 0; |
| 1028 | |
| 1029 | } |
| 1030 | |
| 1031 | |
| 1032 | void Detector::MakeAliceCurrent(Int_t AlignResiduals, Bool_t flagTPC) { |
| 1033 | |
| 1034 | // Numbers taken from |
| 1035 | // 2010 JINST 5 P03003 - Alignment of the ALICE Inner Tracking System with cosmic-ray tracks |
| 1036 | // number for misalingment: private communication with Andrea Dainese |
| 1037 | |
| 1038 | AddLayer((char*)"bpipe",2.94,0.0022); // beam pipe |
| 1039 | AddLayer((char*)"vertex", 0, 0); // dummy vertex for matrix calculation |
| 1040 | AddLayer((char*)"tshld1",11.5,0.0065); // Thermal shield // 1.3% /2 |
| 1041 | AddLayer((char*)"tshld2",31.0,0.0065); // Thermal shield // 1.3% /2 |
| 1042 | AddLayer((char*)"IFC", 77.8,0.01367); // Inner Field cage |
| 1043 | |
| 1044 | if (flagTPC) AddTPC(0.1,0.1); // TPC |
| 1045 | |
| 1046 | // Adding the ITS - current configuration |
| 1047 | |
| 1048 | if (AlignResiduals==0) { |
| 1049 | |
| 1050 | AddLayer((char*)"spd1", 3.9, 0.0114, 0.0012, 0.0100); |
| 1051 | AddLayer((char*)"spd2", 7.6, 0.0114, 0.0012, 0.0100); |
| 1052 | AddLayer((char*)"sdd1",15.0, 0.0113, 0.0035, 0.0025); |
| 1053 | AddLayer((char*)"sdd2",23.9, 0.0126, 0.0035, 0.0025); |
| 1054 | AddLayer((char*)"ssd1",38.0, 0.0083, 0.0020, 0.0830); |
| 1055 | AddLayer((char*)"ssd2",43.0, 0.0086, 0.0020, 0.0830); |
| 1056 | |
| 1057 | } else if (AlignResiduals==1) { |
| 1058 | |
| 1059 | // tracking errors ... |
| 1060 | // (Additional systematic errors due to misalignments) ... |
| 1061 | // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm] |
| 1062 | // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000); |
| 1063 | |
| 1064 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010), |
| 1065 | TMath::Sqrt(0.0100*0.0100+0.0050*0.0050)); |
| 1066 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030), |
| 1067 | TMath::Sqrt(0.0100*0.0100+0.0050*0.0050)); |
| 1068 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500), |
| 1069 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); |
| 1070 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500), |
| 1071 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); |
| 1072 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), |
| 1073 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); |
| 1074 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020), |
| 1075 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); |
| 1076 | |
| 1077 | } else if (AlignResiduals==2) { |
| 1078 | |
| 1079 | // tracking errors ... PLUS ... module misalignment |
| 1080 | |
| 1081 | // itsRecoParam->SetClusterMisalErrorYBOn(0.0010,0.0030,0.0500,0.0500,0.0020,0.0020); // [cm] |
| 1082 | // itsRecoParam->SetClusterMisalErrorZBOn(0.0050,0.0050,0.0050,0.0050,0.1000,0.1000); |
| 1083 | |
| 1084 | // the ITS modules are misalignment with small gaussian smearings with |
| 1085 | // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD |
| 1086 | |
| 1087 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0010*0.0010+0.0008*0.0008), |
| 1088 | TMath::Sqrt(0.0100*0.0100+0.0050*0.0050)); |
| 1089 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0030*0.0030+0.0008*0.0008), |
| 1090 | TMath::Sqrt(0.0100*0.0100+0.0050*0.0050)); |
| 1091 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010), |
| 1092 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); |
| 1093 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0500*0.0500+0.0010*0.0010), |
| 1094 | TMath::Sqrt(0.0025*0.0025+0.0050*0.0050)); |
| 1095 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010), |
| 1096 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); |
| 1097 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0020*0.0020+0.0010*0.0010), |
| 1098 | TMath::Sqrt(0.0830*0.0830+0.1000*0.1000)); |
| 1099 | |
| 1100 | } else { |
| 1101 | |
| 1102 | // the ITS modules are misalignment with small gaussian smearings with |
| 1103 | // sigmarphi ~ 8, 10, 10 micron in SPD, SDD, SSD |
| 1104 | // unknown in Z ???? |
| 1105 | |
| 1106 | AddLayer((char*)"spd1", 3.9, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), |
| 1107 | TMath::Sqrt(0.0100*0.0100+0.000*0.000)); |
| 1108 | AddLayer((char*)"spd2", 7.6, 0.0114, TMath::Sqrt(0.0012*0.0012+0.0008*0.0008), |
| 1109 | TMath::Sqrt(0.0100*0.0100+0.000*0.000)); |
| 1110 | AddLayer((char*)"sdd1",15.0, 0.0113, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010), |
| 1111 | TMath::Sqrt(0.0025*0.0025+0.000*0.000)); |
| 1112 | AddLayer((char*)"sdd2",23.9, 0.0126, TMath::Sqrt(0.0035*0.0035+0.0010*0.0010), |
| 1113 | TMath::Sqrt(0.0025*0.0025+0.000*0.000)); |
| 1114 | AddLayer((char*)"ssd1",38.0, 0.0083, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), |
| 1115 | TMath::Sqrt(0.0830*0.0830+0.000*0.000)); |
| 1116 | AddLayer((char*)"ssd2",43.0, 0.0086, TMath::Sqrt(0.0020*0.0020+0.0010*0.0010), |
| 1117 | TMath::Sqrt(0.0830*0.0830+0.000*0.000)); |
| 1118 | |
| 1119 | |
| 1120 | } |
| 1121 | |
| 1122 | } |
git://git.uio.no / u / mrichter / AliRoot.git / blame