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