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