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