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