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