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 }