]>
Commit | Line | Data |
---|---|---|
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 | } |