Do not include from subdirectories
[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTriggerEff.cxx
CommitLineData
6255180c 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
a42548b0 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().
e0003ae3 26// Author: Pietro Cortese (Universita' del Piemonte Orientale - Alessandria
27// and INFN of Torino)
6255180c 28
490243f3 29
30#include "AliFastMuonTriggerEff.h"
2724a683 31#include "TROOT.h"
32#include "TFile.h"
33#include "stdlib.h"
34#include "TH3.h"
35#include "TObjString.h"
36
490243f3 37#define PLIN printf("%s: %d: ",__FILE__,__LINE__)
6255180c 38
ecda8707 39ClassImp(AliFastMuonTriggerEff)
6255180c 40
41AliFastMuonTriggerEff::AliFastMuonTriggerEff():
e6e76983 42 AliFastResponse("Efficiency", "Muon Trigger Efficiency"),
43 fPtMin(0.),
44 fPtMax(0.),
45 fDpt(0.),
46 fnptb(0),
47 fPhiMin(0.),
48 fPhiMax(0.),
49 fDphi(0.),
50 fnphib(0),
51 fThetaMin(0.),
52 fThetaMax(0.),
53 fDtheta(0.),
54 fnthetab(0),
55 fCut(kLow),
56 fZones(0),
57 fhEffAPt(0),
58 fhEffLPt(0),
59 fhEffHPt(0),
60 fhLX(0),
61 fhLY(0),
62 fhLZ(0),
63 fBkg(0.),
64 fTableTitle(0),
65 fDescription(0),
66 fInt(0),
67 fibx(0),
68 fiby(0),
69 fibz(0)
6255180c 70{
71//
72// Default constructor
73//
e6e76983 74}
75
76AliFastMuonTriggerEff::AliFastMuonTriggerEff(const char* Name, const char* Title):
77 AliFastResponse(Name, Title),
78 fPtMin(0.),
79 fPtMax(0.),
80 fDpt(0.),
81 fnptb(0),
82 fPhiMin(0.),
83 fPhiMax(0.),
84 fDphi(0.),
85 fnphib(0),
86 fThetaMin(0.),
87 fThetaMax(0.),
88 fDtheta(0.),
89 fnthetab(0),
90 fCut(kLow),
91 fZones(0),
92 fhEffAPt(0),
93 fhEffLPt(0),
94 fhEffHPt(0),
95 fhLX(0),
96 fhLY(0),
97 fhLZ(0),
98 fBkg(0.),
99 fTableTitle(0),
100 fDescription(0),
101 fInt(0),
102 fibx(0),
103 fiby(0),
104 fibz(0)
105{
106// Another constructor
6255180c 107}
108
a42548b0 109AliFastMuonTriggerEff::AliFastMuonTriggerEff(const AliFastMuonTriggerEff& eff)
e6e76983 110 :AliFastResponse(eff),
111 fPtMin(0.),
112 fPtMax(0.),
113 fDpt(0.),
114 fnptb(0),
115 fPhiMin(0.),
116 fPhiMax(0.),
117 fDphi(0.),
118 fnphib(0),
119 fThetaMin(0.),
120 fThetaMax(0.),
121 fDtheta(0.),
122 fnthetab(0),
123 fCut(kLow),
124 fZones(0),
125 fhEffAPt(0),
126 fhEffLPt(0),
127 fhEffHPt(0),
128 fhLX(0),
129 fhLY(0),
130 fhLZ(0),
131 fBkg(0.),
132 fTableTitle(0),
133 fDescription(0),
134 fInt(0),
135 fibx(0),
136 fiby(0),
137 fibz(0)
a42548b0 138{
139// Copy constructor
140 eff.Copy(*this);
141}
142
490243f3 143void AliFastMuonTriggerEff::SetCut(Int_t cut)
144{
2724a683 145 //
146 // Set the pt cut
490243f3 147 if(cut==kLow){
148 printf("Selecting Low Pt cut\n");
149 }else if(cut==kHigh){
150 printf("Selecting High Pt cut\n");
151 }else if(cut==kAny){
152 printf("Selecting Lowest Pt cut\n");
153 }else{
154 printf("Don't know cut %d! Selecting Low Pt cut\n",cut);
155 cut=kLow;
156 }
157 fCut = cut;
158
159}
6255180c 160
490243f3 161Int_t AliFastMuonTriggerEff::SetBkgLevel(Float_t Bkg)
6255180c 162{
2724a683 163 //
164 // Set the background level
165 //
490243f3 166 if((Bkg!=0.)) {
167 printf("%s: Warning: requested Bkg: %f\n",
168 __FILE__,Bkg);
169 fBkg=0.;
170 printf("A consistent treatement of the trigger probability\n");
171 printf("within the framework of the fast simulation requires\n");
172 printf("requires background 0\n");
173 printf("%s: fBkg: set to %f\n",
174 __FILE__,fBkg);
175 } else {
176 fBkg=Bkg;
177 }
178 if(fZones!=0.) {
179 Init();
180 }
181 return 0;
182}
6255180c 183
490243f3 184Int_t AliFastMuonTriggerEff::ForceBkgLevel(Float_t Bkg)
185{
2724a683 186 //
187 // Check and enforce consistency of the background level
188 //
490243f3 189 if((Bkg!=0.)) {
190 printf("%s: Warning: requested Bkg: %f\n",
191 __FILE__,Bkg);
192 printf("A consistent treatement of the trigger probability\n");
193 printf("within the framework of the fast simulation\n");
194 printf("requires background 0");
195 printf("%s: Continue with fBkg: %f\n",
196 __FILE__,Bkg);
197 }
198 fBkg=Bkg;
199 if(fZones!=0.) {
200 Init();
201 }
202 return 0;
203}
6255180c 204
c4c8e047 205Int_t AliFastMuonTriggerEff::LoadTables(const Char_t *namet){
2724a683 206 //
207 // Load the trigger tables
208 //
490243f3 209 Char_t hNameA[100],hNameL[100],hNameH[100];
bae6da3c 210 snprintf(hNameA, 100, "hEffAPt%s",namet);
211 snprintf(hNameL, 100, "hEffLPt%s",namet);
212 snprintf(hNameH, 100, "hEffHPt%s",namet);
490243f3 213 fhEffAPt = (TH3F*)gDirectory->Get(hNameA);
214 fhEffLPt = (TH3F*)gDirectory->Get(hNameL);
215 fhEffHPt = (TH3F*)gDirectory->Get(hNameH);
216 if(!fhEffAPt){
217 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameA);
218 return -1;
219 }
220 if(!fhEffLPt){
221 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameL);
222 return -2;
223 }
224 if(!fhEffHPt){
225 PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameH);
226 return -3;
227 }
228 return 0;
6255180c 229}
230
490243f3 231void AliFastMuonTriggerEff::Init()
6255180c 232{
233//
490243f3 234// Initialization
6255180c 235//
490243f3 236 fZones=0;
237 Char_t file[100]="$(ALICE_ROOT)/FASTSIM/data/MUONtriggerLUT_V2.4nvdn.root";
238 printf("Initializing %s / %s\n", fName.Data(), fTitle.Data());
239 printf("using data from file: %s\n",file);
240 printf("AliFastMuonTriggerEff: Initialization with background level: %f\n",fBkg);
6255180c 241 TFile *f = new TFile(file);
242 if(f->IsZombie()) {
490243f3 243 PLIN; printf("Cannot open file: %s\n",file);
6255180c 244 return;
245 }
246 f->ls();
490243f3 247 Int_t intb=0;
248 Char_t namet[10];
249 if(TMath::Abs(fBkg)<0.00001){
bae6da3c 250 snprintf(namet, 10, "00");
490243f3 251 }else if(TMath::Abs(fBkg-0.5)<0.00001){
bae6da3c 252 snprintf(namet, 10, "05");
490243f3 253 }else if(TMath::Abs(fBkg-1.0)<0.00001){
bae6da3c 254 snprintf(namet, 10, "10");
490243f3 255 }else{
256 PLIN; printf("A table for Bkg level: %f does not exists\n",fBkg);
257 intb=1;
6255180c 258 }
490243f3 259 if(intb){ // Interpolation between background levels
260 PLIN; printf("Interpolating Bkg level: %f\n",fBkg);
261 TH3F* ha1,*hl1,*hh1,*ha2,*hl2,*hh2,*ha0,*hl0,*hh0;
262 Char_t name1[10],name2[10]; Float_t b1,b2;
263 if(fBkg>0&&fBkg<0.5){
bae6da3c 264 snprintf(name1,10, "00");
265 snprintf(name2,10, "05");
490243f3 266 b1=0.;
267 b2=0.5;
268 }else if(fBkg>0.5){
bae6da3c 269 snprintf(name1, 10, "05");
270 snprintf(name2, 10, "10");
490243f3 271 b1=0.5;
272 b2=1.0;
273 if(fBkg>1.0){
274 for(Int_t i=0; i<10;i++){
275 PLIN; printf("WARNING!!!! You are extrapolating above background 1.0\n");
276 }
6255180c 277 }
490243f3 278 }else{
279 PLIN; printf("Bkg level: %f is not supported\n",fBkg);
280 return;
281 }
282 if(LoadTables(name1)){
283 PLIN; printf("Error in loading trigger tables\n");
284 return;
285 }
286 PLIN; printf("We use tables for %f and %f to interpolate %f Bkg level\n",b1,b2,fBkg);
287 ha0=(TH3F*)fhEffAPt->Clone("hEffAPtXX"); ha0->Reset();
288 hl0=(TH3F*)fhEffLPt->Clone("hEffLPtXX"); hl0->Reset();
289 hh0=(TH3F*)fhEffHPt->Clone("hEffHPtXX"); hh0->Reset();
290 ha1=fhEffAPt;
291 hl1=fhEffLPt;
292 hh1=fhEffHPt;
293 if(LoadTables(name2)){
294 PLIN; printf("Error in loading trigger tables\n");
295 return;
296 }
297 ha2=fhEffAPt;
298 hl2=fhEffLPt;
299 hh2=fhEffHPt;
300 fhEffAPt=ha0;
301 fhEffLPt=hl0;
302 fhEffHPt=hh0;
303 Int_t nnx=ha0->GetNbinsX()+1;
304 Int_t nny=ha0->GetNbinsY()+1;
305 Int_t nnz=ha0->GetNbinsZ()+1;
306 for(Int_t ix=0; ix<=nnx; ix++){
307 for(Int_t iy=0; iy<=nny; iy++){
308 for(Int_t iz=0; iz<=nnz; iz++){
309 Double_t y1,y2; Float_t cont;
310 y1=ha1->GetBinContent(ix,iy,iz); y2=ha2->GetBinContent(ix,iy,iz);
311 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
312 fhEffAPt->SetBinContent(ix,iy,iz,cont);
313 y1=hl1->GetBinContent(ix,iy,iz); y2=hl2->GetBinContent(ix,iy,iz);
314 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
315 fhEffLPt->SetBinContent(ix,iy,iz,cont);
316 y1=hh1->GetBinContent(ix,iy,iz); y2=hh2->GetBinContent(ix,iy,iz);
317 cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
318 fhEffHPt->SetBinContent(ix,iy,iz,cont);
319 }
6255180c 320 }
490243f3 321 }
322 }else{ // Use tables computed for selected backgound levels
323 printf("Loading tables for background level: %f\n",fBkg);
324 if(LoadTables(namet)){
325 PLIN; printf("Error in loading trigger tables\n");
326 return;
327 }
6255180c 328 }
490243f3 329 fhEffAPt->SetDirectory(0);
330 fhEffLPt->SetDirectory(0);
331 fhEffHPt->SetDirectory(0);
332 fhLX=fhEffLPt->GetXaxis();
333 fhLY=fhEffLPt->GetYaxis();
334 fhLZ=fhEffLPt->GetZaxis();
335//
336//
337 if(f->Get("Description"))
338 {
339 fDescription=((TObjString*)f->Get("Description"))->GetString();
340 printf("%s\n",fDescription.Data());
6255180c 341 }
490243f3 342
343 fThetaMin = fhEffLPt->GetXaxis()->GetXmin();
344 fThetaMax = fhEffLPt->GetXaxis()->GetXmax();
345 fnthetab=fhEffLPt->GetNbinsX();
346 fDtheta = (fThetaMax-fThetaMin)/fnthetab;
347
348 fPhiMin = fhEffLPt->GetYaxis()->GetXmin();
349 fPhiMax = fhEffLPt->GetYaxis()->GetXmax();
350 fnphib=fhEffLPt->GetNbinsY();
351 fDphi = (fPhiMax-fPhiMin)/fnphib;
352
353 fPtMin=fhEffLPt->GetZaxis()->GetXmin();
354 fPtMax=fhEffLPt->GetZaxis()->GetXmax();
355 fnptb=fhEffLPt->GetNbinsZ();
356 fDpt = (fPtMax-fPtMin)/fnptb;
357
e0003ae3 358 printf("***** This version of AliFastMuonTriggerEff can use both *****\n");
359 printf("***** new and old ALICE reference frames depending on *****\n");
360 printf("***** which LUT has been loaded. You can find below some *****\n");
361 printf("***** information on the current parametrization: *****\n");
490243f3 362 printf("%4d bins in theta [%f:%f]\n",fnthetab,fThetaMin,fThetaMax);
363 printf("%4d bins in phi [%f:%f]\n",fnphib,fPhiMin,fPhiMax);
364 printf("%4d bins in pt [%f:%f]\n",fnptb,fPtMin,fPtMax);
365
366 fZones=fnthetab*fnphib;
367
6255180c 368 f->Close();
490243f3 369 if(fInt==0) {
370 printf("Interpolation of trigger efficiencies is off!\n");
371 } else {
372 printf("Interpolation of trigger efficiencies is on!\n");
373 }
6255180c 374}
375
490243f3 376void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,Float_t theta,
377 Float_t phi, Float_t& effLow, Float_t& effHigh, Float_t& effAny)
6255180c 378{
490243f3 379 //
380 // Trigger efficiency for pt, theta, phi (low, high and "any" cut)
381 //
382#ifdef MYTRIGDEBUG
383 printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
384#endif
6255180c 385 effLow=0.;
386 effHigh=0.;
490243f3 387 effAny=0;
6255180c 388 if(fZones==0) {
389 printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
390 return;
391 }
392 if(pt<0) {
393 printf("Warning: pt: %f < 0. GeV/c\n",pt);
394 return;
395 }
490243f3 396
397 Int_t iPt = fhLZ->FindBin((Double_t)pt);
398 if(iPt>fnptb)iPt=fnptb;
6255180c 399 Int_t iPhi = Int_t((phi-fPhiMin)/fDphi);
400 if(phi<fPhiMin)iPhi=iPhi-1;
490243f3 401 Int_t iTheta = fhLX->FindBin((Double_t)theta);
6255180c 402#ifdef MYTRIGDEBUG
490243f3 403 printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
404 printf(" 0:%1d iPt iTheta iPhi: %d %d %d\n",fInt,iPt,iTheta,iPhi);
6255180c 405#endif
490243f3 406 iPhi=iPhi-2*fnphib*(iPhi/(2*fnphib));
6255180c 407#ifdef MYTRIGDEBUG
490243f3 408 printf(" 1:%1d iPhi converted to: %d for angle equivalence\n",fInt,iPhi);
6255180c 409#endif
410 if(iPhi<0)iPhi=-iPhi-1;
490243f3 411 if(iPhi>(fnphib-1))iPhi=2*fnphib-1-iPhi;
6255180c 412#ifdef MYTRIGDEBUG
490243f3 413 printf(" 2:%1d iPhi converted to: %d for the symmetry of the spectrometer\n",fInt,iPhi);
6255180c 414#endif
415 if(charge==1.){
416 } else if(charge==-1.) {
490243f3 417 iPhi=fnphib-1-iPhi;
6255180c 418#ifdef MYTRIGDEBUG
490243f3 419 printf(" 3:%1d iPhi converted to: %d for the charge symmetry\n",fInt,iPhi);
6255180c 420#endif
421 } else {
422 printf("Warning: not understand charge: %f\n",charge);
423 return;
424 }
490243f3 425 if(iTheta<=0||iTheta>fnthetab) {
6255180c 426 printf("Warning: theta: %f outside acceptance\n",theta);
427 return;
428 }
429 if(iPt<0) {
430 printf("Warning: what do you mean with pt: %f <0?\n",pt);
431 return;
432 }
490243f3 433 iPhi++;
6255180c 434#ifdef MYTRIGDEBUG
490243f3 435 printf(" 4:%1d Getting: iTheta, iPhi, iPt: %d %d %d\n",
436 fInt,iTheta,iPhi,iPt);
6255180c 437#endif
490243f3 438 effLow =fhEffLPt->GetBinContent(iTheta,iPhi,iPt);
439 effHigh=fhEffHPt->GetBinContent(iTheta,iPhi,iPt);
440 effAny =fhEffAPt->GetBinContent(iTheta,iPhi,iPt);
6255180c 441#ifdef MYTRIGDEBUG
490243f3 442 printf(" 4:%1d Result: charge, iTheta, iPhi, iPt: %f %d %d %d effLow: %f, effHigh: %f, effAny: %f\n",
443 fInt,charge,iTheta,iPhi,iPt,effLow,effHigh,effAny);
6255180c 444#endif
490243f3 445
446 if(fInt==1) {
447 Float_t angl,angh,anga;
448 Float_t effLowp,effHighp,effAnyp;
449 Float_t ptc=(iPt+0.5)*fDpt; // The center of current bin
450 #ifdef MYTRIGDEBUG
451 printf(" 5:1 The center of current bin iPt: %d is: %f\n",iPt,ptc);
452 #endif
453 if(iPt==fnptb) {
454 #ifdef MYTRIGDEBUG
455 printf(" 6:1 No more points above! No interpolation is needed!\n");
456 #endif
457 return;
620294e8 458 }else if(TMath::Abs(ptc-pt) < 1.e-10){
490243f3 459 #ifdef MYTRIGDEBUG
460 printf(" 6:1 No interpolation is needed!\n");
461 #endif
462 return;
463 }else if(ptc>pt){
464 // Looking for previous point
465 if(iPt>1) {
466 effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt-1);
467 effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt-1);
468 effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt-1);
469 #ifdef MYTRIGDEBUG
470 printf(" 7:1 A simple look to previous point: %d: %f %f\n",iPt-1,effLowp,effHighp);
471 #endif
472 } else {
473 effLowp=0.;
474 effHighp=0.;
475 effAnyp=0;
476 #ifdef MYTRIGDEBUG
477 printf(" 8:1 result is: %f %f %f\n",effLowp,effHighp,effAnyp);
478 #endif
479 }
480 angl=(effLow-effLowp)/fDpt;
481 angh=(effHigh-effHighp)/fDpt;
482 anga=(effAny-effAnyp)/fDpt;
483 }else{
484 // Looking for next point
485 if(iPt<fnptb) {
486 effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt+1);
487 effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt+1);
488 effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt+1);
489 #ifdef MYTRIGDEBUG
490 printf(" 7:1 A simple look to next point: %d: %f %f %f\n",iPt-1,effLowp,effHighp,effAnyp);
491 #endif
492 } else {
493 effLowp=effLow;
494 effHighp=effHigh;
495 effAnyp=effAny;
496 #ifdef MYTRIGDEBUG
497 printf(" 8:1 result is: pt: %f %f %f\n",effLowp,effHighp,effAnyp);
498 #endif
499 }
500 angl=(effLowp-effLow)/fDpt;
501 angh=(effHighp-effHigh)/fDpt;
502 anga=(effAnyp-effAny)/fDpt;
503 }
504 effLow=effLow+angl*(pt-ptc);
505 effHigh=effHigh+angh*(pt-ptc);
506 effAny=effAny+anga*(pt-ptc);
507 #ifdef MYTRIGDEBUG
508 printf(" 9:1 the interpolation coefficients are: %f %f %f\n",angl,angh,anga);
509 #endif
510 }
511 #ifdef MYTRIGDEBUG
512 printf("10:%1d effLow, effHigh=%f %f %f\n",fInt,effLow,effHigh,effAny);
513 #endif
514 return;
515}
516
bc77d57e 517
6255180c 518
519Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
520 Float_t theta, Float_t phi)
521{
522 //
523 // Trigger efficiency for pt, theta, phi depending of fCut
524 //
525 if(fZones==0) {
526 printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
527 return 0.;
528 }
529 Float_t eff;
490243f3 530 Float_t effLow, effHigh, effAny;
6255180c 531
490243f3 532 Evaluate(charge,pt,theta,phi,effLow,effHigh,effAny);
6255180c 533 if (fCut == kLow)
534 eff = effLow;
490243f3 535 else if (fCut == kHigh)
6255180c 536 eff = effHigh;
490243f3 537 else if (fCut == kAny)
538 eff = effAny;
539 else
540 eff = 0;
6255180c 541
542 return eff;
543}
a42548b0 544
545AliFastMuonTriggerEff& AliFastMuonTriggerEff::operator=(const AliFastMuonTriggerEff& rhs)
546{
547// Assignment operator
548 rhs.Copy(*this);
549 return *this;
550}
551