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