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