Geometry builder classes moved from base to sim.
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTriggerEff.cxx
... / ...
CommitLineData
1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
17// The trigger parametrization is computed for background levels 0., 0.5 and 1.
18// In order to set a background level different from 0 it is necessary to
19// explicitly force it with:
20// ForceBkgLevel(BkgLevel).
21// For intermediate background levels, the trigger response is linearly
22// interpolated between these values.
23// There is increased granularity in the pT region below 3 GeV. Although
24// it does not seem to be necessary it is also possible to interpolate
25// between pT bins using SetInt().
26
27
28#include "AliFastMuonTriggerEff.h"
29#include "TROOT.h"
30#include "TFile.h"
31#include "stdlib.h"
32#include "TH3.h"
33#include "TObjString.h"
34
35#define PLIN printf("%s: %d: ",__FILE__,__LINE__)
36
37ClassImp(AliFastMuonTriggerEff)
38
39AliFastMuonTriggerEff::AliFastMuonTriggerEff():
40 AliFastResponse("Efficiency", "Muon Trigger Efficiency")
41{
42//
43// Default constructor
44//
45 SetCut(kLow);
46 fZones=0;
47 SetBkgLevel(0.);
48 UnsetInt();
49}
50
51AliFastMuonTriggerEff::AliFastMuonTriggerEff(const AliFastMuonTriggerEff& eff)
52 :AliFastResponse(eff)
53{
54// Copy constructor
55 eff.Copy(*this);
56}
57
58void AliFastMuonTriggerEff::SetCut(Int_t cut)
59{
60 //
61 // Set the pt cut
62 if(cut==kLow){
63 printf("Selecting Low Pt cut\n");
64 }else if(cut==kHigh){
65 printf("Selecting High Pt cut\n");
66 }else if(cut==kAny){
67 printf("Selecting Lowest Pt cut\n");
68 }else{
69 printf("Don't know cut %d! Selecting Low Pt cut\n",cut);
70 cut=kLow;
71 }
72 fCut = cut;
73
74}
75
76Int_t AliFastMuonTriggerEff::SetBkgLevel(Float_t Bkg)
77{
78 //
79 // Set the background level
80 //
81 if((Bkg!=0.)) {
82 printf("%s: Warning: requested Bkg: %f\n",
83 __FILE__,Bkg);
84 fBkg=0.;
85 printf("A consistent treatement of the trigger probability\n");
86 printf("within the framework of the fast simulation requires\n");
87 printf("requires background 0\n");
88 printf("%s: fBkg: set to %f\n",
89 __FILE__,fBkg);
90 } else {
91 fBkg=Bkg;
92 }
93 if(fZones!=0.) {
94 Init();
95 }
96 return 0;
97}
98
99Int_t AliFastMuonTriggerEff::ForceBkgLevel(Float_t Bkg)
100{
101 //
102 // Check and enforce consistency of the background level
103 //
104 if((Bkg!=0.)) {
105 printf("%s: Warning: requested Bkg: %f\n",
106 __FILE__,Bkg);
107 printf("A consistent treatement of the trigger probability\n");
108 printf("within the framework of the fast simulation\n");
109 printf("requires background 0");
110 printf("%s: Continue with fBkg: %f\n",
111 __FILE__,Bkg);
112 }
113 fBkg=Bkg;
114 if(fZones!=0.) {
115 Init();
116 }
117 return 0;
118}
119
120Int_t AliFastMuonTriggerEff::LoadTables(Char_t *namet=""){
121 //
122 // Load the trigger tables
123 //
124 Char_t hNameA[100],hNameL[100],hNameH[100];
125 sprintf(hNameA,"hEffAPt%s",namet);
126 sprintf(hNameL,"hEffLPt%s",namet);
127 sprintf(hNameH,"hEffHPt%s",namet);
128 fhEffAPt = (TH3F*)gDirectory->Get(hNameA);
129 fhEffLPt = (TH3F*)gDirectory->Get(hNameL);
130 fhEffHPt = (TH3F*)gDirectory->Get(hNameH);
131 if(!fhEffAPt){
132 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameA);
133 return -1;
134 }
135 if(!fhEffLPt){
136 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameL);
137 return -2;
138 }
139 if(!fhEffHPt){
140 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameH);
141 return -3;
142 }
143 return 0;
144}
145
146void AliFastMuonTriggerEff::Init()
147{
148//
149// Initialization
150//
151 fZones=0;
152 Char_t file[100]="$(ALICE_ROOT)/FASTSIM/data/MUONtriggerLUT_V2.4nvdn.root";
153 printf("Initializing %s / %s\n", fName.Data(), fTitle.Data());
154 printf("using data from file: %s\n",file);
155 printf("AliFastMuonTriggerEff: Initialization with background level: %f\n",fBkg);
156 TFile *f = new TFile(file);
157 if(f->IsZombie()) {
158 PLIN; printf("Cannot open file: %s\n",file);
159 return;
160 }
161 f->ls();
162 Int_t intb=0;
163 Char_t namet[10];
164 if(TMath::Abs(fBkg)<0.00001){
165 sprintf(namet,"00");
166 }else if(TMath::Abs(fBkg-0.5)<0.00001){
167 sprintf(namet,"05");
168 }else if(TMath::Abs(fBkg-1.0)<0.00001){
169 sprintf(namet,"10");
170 }else{
171 PLIN; printf("A table for Bkg level: %f does not exists\n",fBkg);
172 intb=1;
173 }
174 if(intb){ // Interpolation between background levels
175 PLIN; printf("Interpolating Bkg level: %f\n",fBkg);
176 TH3F* ha1,*hl1,*hh1,*ha2,*hl2,*hh2,*ha0,*hl0,*hh0;
177 Char_t name1[10],name2[10]; Float_t b1,b2;
178 if(fBkg>0&&fBkg<0.5){
179 sprintf(name1,"00");
180 sprintf(name2,"05");
181 b1=0.;
182 b2=0.5;
183 }else if(fBkg>0.5){
184 sprintf(name1,"05");
185 sprintf(name2,"10");
186 b1=0.5;
187 b2=1.0;
188 if(fBkg>1.0){
189 for(Int_t i=0; i<10;i++){
190 PLIN; printf("WARNING!!!! You are extrapolating above background 1.0\n");
191 }
192 }
193 }else{
194 PLIN; printf("Bkg level: %f is not supported\n",fBkg);
195 return;
196 }
197 if(LoadTables(name1)){
198 PLIN; printf("Error in loading trigger tables\n");
199 return;
200 }
201 PLIN; printf("We use tables for %f and %f to interpolate %f Bkg level\n",b1,b2,fBkg);
202 ha0=(TH3F*)fhEffAPt->Clone("hEffAPtXX"); ha0->Reset();
203 hl0=(TH3F*)fhEffLPt->Clone("hEffLPtXX"); hl0->Reset();
204 hh0=(TH3F*)fhEffHPt->Clone("hEffHPtXX"); hh0->Reset();
205 ha1=fhEffAPt;
206 hl1=fhEffLPt;
207 hh1=fhEffHPt;
208 if(LoadTables(name2)){
209 PLIN; printf("Error in loading trigger tables\n");
210 return;
211 }
212 ha2=fhEffAPt;
213 hl2=fhEffLPt;
214 hh2=fhEffHPt;
215 fhEffAPt=ha0;
216 fhEffLPt=hl0;
217 fhEffHPt=hh0;
218 Int_t nnx=ha0->GetNbinsX()+1;
219 Int_t nny=ha0->GetNbinsY()+1;
220 Int_t nnz=ha0->GetNbinsZ()+1;
221 for(Int_t ix=0; ix<=nnx; ix++){
222 for(Int_t iy=0; iy<=nny; iy++){
223 for(Int_t iz=0; iz<=nnz; iz++){
224 Double_t y1,y2; Float_t cont;
225 y1=ha1->GetBinContent(ix,iy,iz); y2=ha2->GetBinContent(ix,iy,iz);
226 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
227 fhEffAPt->SetBinContent(ix,iy,iz,cont);
228 y1=hl1->GetBinContent(ix,iy,iz); y2=hl2->GetBinContent(ix,iy,iz);
229 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
230 fhEffLPt->SetBinContent(ix,iy,iz,cont);
231 y1=hh1->GetBinContent(ix,iy,iz); y2=hh2->GetBinContent(ix,iy,iz);
232 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
233 fhEffHPt->SetBinContent(ix,iy,iz,cont);
234 }
235 }
236 }
237 }else{ // Use tables computed for selected backgound levels
238 printf("Loading tables for background level: %f\n",fBkg);
239 if(LoadTables(namet)){
240 PLIN; printf("Error in loading trigger tables\n");
241 return;
242 }
243 }
244 fhEffAPt->SetDirectory(0);
245 fhEffLPt->SetDirectory(0);
246 fhEffHPt->SetDirectory(0);
247 fhLX=fhEffLPt->GetXaxis();
248 fhLY=fhEffLPt->GetYaxis();
249 fhLZ=fhEffLPt->GetZaxis();
250//
251//
252 if(f->Get("Description"))
253 {
254 fDescription=((TObjString*)f->Get("Description"))->GetString();
255 printf("%s\n",fDescription.Data());
256 }
257
258 fThetaMin = fhEffLPt->GetXaxis()->GetXmin();
259 fThetaMax = fhEffLPt->GetXaxis()->GetXmax();
260 fnthetab=fhEffLPt->GetNbinsX();
261 fDtheta = (fThetaMax-fThetaMin)/fnthetab;
262
263 fPhiMin = fhEffLPt->GetYaxis()->GetXmin();
264 fPhiMax = fhEffLPt->GetYaxis()->GetXmax();
265 fnphib=fhEffLPt->GetNbinsY();
266 fDphi = (fPhiMax-fPhiMin)/fnphib;
267
268 fPtMin=fhEffLPt->GetZaxis()->GetXmin();
269 fPtMax=fhEffLPt->GetZaxis()->GetXmax();
270 fnptb=fhEffLPt->GetNbinsZ();
271 fDpt = (fPtMax-fPtMin)/fnptb;
272
273 printf("%4d bins in theta [%f:%f]\n",fnthetab,fThetaMin,fThetaMax);
274 printf("%4d bins in phi [%f:%f]\n",fnphib,fPhiMin,fPhiMax);
275 printf("%4d bins in pt [%f:%f]\n",fnptb,fPtMin,fPtMax);
276
277 fZones=fnthetab*fnphib;
278
279 f->Close();
280 if(fInt==0) {
281 printf("Interpolation of trigger efficiencies is off!\n");
282 } else {
283 printf("Interpolation of trigger efficiencies is on!\n");
284 }
285}
286
287void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,Float_t theta,
288 Float_t phi, Float_t& effLow, Float_t& effHigh, Float_t& effAny)
289{
290 //
291 // Trigger efficiency for pt, theta, phi (low, high and "any" cut)
292 //
293#ifdef MYTRIGDEBUG
294 printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
295#endif
296 effLow=0.;
297 effHigh=0.;
298 effAny=0;
299 if(fZones==0) {
300 printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
301 return;
302 }
303 if(pt<0) {
304 printf("Warning: pt: %f < 0. GeV/c\n",pt);
305 return;
306 }
307
308 Int_t iPt = fhLZ->FindBin((Double_t)pt);
309 if(iPt>fnptb)iPt=fnptb;
310 Int_t iPhi = Int_t((phi-fPhiMin)/fDphi);
311 if(phi<fPhiMin)iPhi=iPhi-1;
312 Int_t iTheta = fhLX->FindBin((Double_t)theta);
313#ifdef MYTRIGDEBUG
314 printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
315 printf(" 0:%1d iPt iTheta iPhi: %d %d %d\n",fInt,iPt,iTheta,iPhi);
316#endif
317 iPhi=iPhi-2*fnphib*(iPhi/(2*fnphib));
318#ifdef MYTRIGDEBUG
319 printf(" 1:%1d iPhi converted to: %d for angle equivalence\n",fInt,iPhi);
320#endif
321 if(iPhi<0)iPhi=-iPhi-1;
322 if(iPhi>(fnphib-1))iPhi=2*fnphib-1-iPhi;
323#ifdef MYTRIGDEBUG
324 printf(" 2:%1d iPhi converted to: %d for the symmetry of the spectrometer\n",fInt,iPhi);
325#endif
326 if(charge==1.){
327 } else if(charge==-1.) {
328 iPhi=fnphib-1-iPhi;
329#ifdef MYTRIGDEBUG
330 printf(" 3:%1d iPhi converted to: %d for the charge symmetry\n",fInt,iPhi);
331#endif
332 } else {
333 printf("Warning: not understand charge: %f\n",charge);
334 return;
335 }
336 if(iTheta<=0||iTheta>fnthetab) {
337 printf("Warning: theta: %f outside acceptance\n",theta);
338 return;
339 }
340 if(iPt<0) {
341 printf("Warning: what do you mean with pt: %f <0?\n",pt);
342 return;
343 }
344 iPhi++;
345#ifdef MYTRIGDEBUG
346 printf(" 4:%1d Getting: iTheta, iPhi, iPt: %d %d %d\n",
347 fInt,iTheta,iPhi,iPt);
348#endif
349 effLow =fhEffLPt->GetBinContent(iTheta,iPhi,iPt);
350 effHigh=fhEffHPt->GetBinContent(iTheta,iPhi,iPt);
351 effAny =fhEffAPt->GetBinContent(iTheta,iPhi,iPt);
352#ifdef MYTRIGDEBUG
353 printf(" 4:%1d Result: charge, iTheta, iPhi, iPt: %f %d %d %d effLow: %f, effHigh: %f, effAny: %f\n",
354 fInt,charge,iTheta,iPhi,iPt,effLow,effHigh,effAny);
355#endif
356
357 if(fInt==1) {
358 Float_t angl,angh,anga;
359 Float_t effLowp,effHighp,effAnyp;
360 Float_t ptc=(iPt+0.5)*fDpt; // The center of current bin
361 #ifdef MYTRIGDEBUG
362 printf(" 5:1 The center of current bin iPt: %d is: %f\n",iPt,ptc);
363 #endif
364 if(iPt==fnptb) {
365 #ifdef MYTRIGDEBUG
366 printf(" 6:1 No more points above! No interpolation is needed!\n");
367 #endif
368 return;
369 }else if(ptc==pt){
370 #ifdef MYTRIGDEBUG
371 printf(" 6:1 No interpolation is needed!\n");
372 #endif
373 return;
374 }else if(ptc>pt){
375 // Looking for previous point
376 if(iPt>1) {
377 effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt-1);
378 effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt-1);
379 effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt-1);
380 #ifdef MYTRIGDEBUG
381 printf(" 7:1 A simple look to previous point: %d: %f %f\n",iPt-1,effLowp,effHighp);
382 #endif
383 } else {
384 effLowp=0.;
385 effHighp=0.;
386 effAnyp=0;
387 #ifdef MYTRIGDEBUG
388 printf(" 8:1 result is: %f %f %f\n",effLowp,effHighp,effAnyp);
389 #endif
390 }
391 angl=(effLow-effLowp)/fDpt;
392 angh=(effHigh-effHighp)/fDpt;
393 anga=(effAny-effAnyp)/fDpt;
394 }else{
395 // Looking for next point
396 if(iPt<fnptb) {
397 effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt+1);
398 effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt+1);
399 effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt+1);
400 #ifdef MYTRIGDEBUG
401 printf(" 7:1 A simple look to next point: %d: %f %f %f\n",iPt-1,effLowp,effHighp,effAnyp);
402 #endif
403 } else {
404 effLowp=effLow;
405 effHighp=effHigh;
406 effAnyp=effAny;
407 #ifdef MYTRIGDEBUG
408 printf(" 8:1 result is: pt: %f %f %f\n",effLowp,effHighp,effAnyp);
409 #endif
410 }
411 angl=(effLowp-effLow)/fDpt;
412 angh=(effHighp-effHigh)/fDpt;
413 anga=(effAnyp-effAny)/fDpt;
414 }
415 effLow=effLow+angl*(pt-ptc);
416 effHigh=effHigh+angh*(pt-ptc);
417 effAny=effAny+anga*(pt-ptc);
418 #ifdef MYTRIGDEBUG
419 printf(" 9:1 the interpolation coefficients are: %f %f %f\n",angl,angh,anga);
420 #endif
421 }
422 #ifdef MYTRIGDEBUG
423 printf("10:%1d effLow, effHigh=%f %f %f\n",fInt,effLow,effHigh,effAny);
424 #endif
425 return;
426}
427
428
429
430Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
431 Float_t theta, Float_t phi)
432{
433 //
434 // Trigger efficiency for pt, theta, phi depending of fCut
435 //
436 if(fZones==0) {
437 printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
438 return 0.;
439 }
440 Float_t eff;
441 Float_t effLow, effHigh, effAny;
442
443 Evaluate(charge,pt,theta,phi,effLow,effHigh,effAny);
444 if (fCut == kLow)
445 eff = effLow;
446 else if (fCut == kHigh)
447 eff = effHigh;
448 else if (fCut == kAny)
449 eff = effAny;
450 else
451 eff = 0;
452
453 return eff;
454}
455
456AliFastMuonTriggerEff& AliFastMuonTriggerEff::operator=(const AliFastMuonTriggerEff& rhs)
457{
458// Assignment operator
459 rhs.Copy(*this);
460 return *this;
461}
462