[u/mrichter/AliRoot.git] / FASTSIM / AliFastMuonTriggerEff.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */
// The trigger parametrization is computed for background levels 0., 0.5 and 1.
// In order to set a background level different from 0 it is necessary to 
// explicitly force it with:
// ForceBkgLevel(BkgLevel).
// For intermediate background levels, the trigger response is linearly 
// interpolated between these values.
// There is increased granularity in the pT region below 3 GeV. Although
// it does not seem to be necessary it is also possible to interpolate
// between pT bins using SetInt().


#include "AliFastMuonTriggerEff.h"
#include "TROOT.h"
#include "TFile.h"
#include "stdlib.h"
#include "TH3.h"
#include "TObjString.h"

#define PLIN printf("%s: %d: ",__FILE__,__LINE__)

ClassImp(AliFastMuonTriggerEff)

AliFastMuonTriggerEff::AliFastMuonTriggerEff():
    AliFastResponse("Efficiency", "Muon Trigger Efficiency"),
    fPtMin(0.),
    fPtMax(0.),
    fDpt(0.),
    fnptb(0),
    fPhiMin(0.),
    fPhiMax(0.),
    fDphi(0.),
    fnphib(0),
    fThetaMin(0.),
    fThetaMax(0.),
    fDtheta(0.),
    fnthetab(0),
    fCut(kLow),
    fZones(0),
    fhEffAPt(0),
    fhEffLPt(0),
    fhEffHPt(0),
    fhLX(0),
    fhLY(0),
    fhLZ(0),
    fBkg(0.),
    fTableTitle(0),
    fDescription(0),
    fInt(0),
    fibx(0),
    fiby(0),
    fibz(0)
{
//
// Default constructor
//
}

AliFastMuonTriggerEff::AliFastMuonTriggerEff(const char* Name, const char* Title):
    AliFastResponse(Name, Title),
    fPtMin(0.),
    fPtMax(0.),
    fDpt(0.),
    fnptb(0),
    fPhiMin(0.),
    fPhiMax(0.),
    fDphi(0.),
    fnphib(0),
    fThetaMin(0.),
    fThetaMax(0.),
    fDtheta(0.),
    fnthetab(0),
    fCut(kLow),
    fZones(0),
    fhEffAPt(0),
    fhEffLPt(0),
    fhEffHPt(0),
    fhLX(0),
    fhLY(0),
    fhLZ(0),
    fBkg(0.),
    fTableTitle(0),
    fDescription(0),
    fInt(0),
    fibx(0),
    fiby(0),
    fibz(0)
{
// Another constructor
}

AliFastMuonTriggerEff::AliFastMuonTriggerEff(const AliFastMuonTriggerEff& eff)
    :AliFastResponse(eff),
    fPtMin(0.),
    fPtMax(0.),
    fDpt(0.),
    fnptb(0),
    fPhiMin(0.),
    fPhiMax(0.),
    fDphi(0.),
    fnphib(0),
    fThetaMin(0.),
    fThetaMax(0.),
    fDtheta(0.),
    fnthetab(0),
    fCut(kLow),
    fZones(0),
    fhEffAPt(0),
    fhEffLPt(0),
    fhEffHPt(0),
    fhLX(0),
    fhLY(0),
    fhLZ(0),
    fBkg(0.),
    fTableTitle(0),
    fDescription(0),
    fInt(0),
    fibx(0),
    fiby(0),
    fibz(0)
{
// Copy constructor
    eff.Copy(*this);
}

void AliFastMuonTriggerEff::SetCut(Int_t cut) 
{  
  //
  // Set the pt cut
  if(cut==kLow){
    printf("Selecting Low Pt cut\n");
  }else if(cut==kHigh){
    printf("Selecting High Pt cut\n");
  }else if(cut==kAny){
    printf("Selecting Lowest Pt cut\n");
  }else{
    printf("Don't know cut %d! Selecting Low Pt cut\n",cut);
    cut=kLow;
  }
  fCut = cut;
  
}

Int_t AliFastMuonTriggerEff::SetBkgLevel(Float_t Bkg)
{
  //
  // Set the background level
  //
    if((Bkg!=0.)) {
       printf("%s: Warning: requested Bkg: %f\n",
       __FILE__,Bkg);
       fBkg=0.;
       printf("A consistent treatement of the trigger probability\n");
       printf("within the framework of the fast simulation requires\n");
       printf("requires background 0\n");
       printf("%s: fBkg: set to %f\n",
       __FILE__,fBkg);
    } else {
      fBkg=Bkg;
    }
    if(fZones!=0.) {
	Init();
    }
    return 0;
}

Int_t AliFastMuonTriggerEff::ForceBkgLevel(Float_t Bkg)
{
  //
  // Check and enforce consistency of the background level 
  // 
    if((Bkg!=0.)) {
       printf("%s: Warning: requested Bkg: %f\n",
       __FILE__,Bkg);
       printf("A consistent treatement of the trigger probability\n");
       printf("within the framework of the fast simulation\n");
       printf("requires background 0");
       printf("%s: Continue with fBkg: %f\n",
       __FILE__,Bkg);
    }
    fBkg=Bkg;
    if(fZones!=0.) {
	Init();
    }
    return 0;
}

Int_t AliFastMuonTriggerEff::LoadTables(Char_t *namet=""){
  //
  // Load the trigger tables
  //
    Char_t hNameA[100],hNameL[100],hNameH[100];
    sprintf(hNameA,"hEffAPt%s",namet);
    sprintf(hNameL,"hEffLPt%s",namet);
    sprintf(hNameH,"hEffHPt%s",namet);
    fhEffAPt = (TH3F*)gDirectory->Get(hNameA);
    fhEffLPt = (TH3F*)gDirectory->Get(hNameL);
    fhEffHPt = (TH3F*)gDirectory->Get(hNameH);
    if(!fhEffAPt){
      PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameA);
      return -1;
    }
    if(!fhEffLPt){
      PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameL);
      return -2;
    }
    if(!fhEffHPt){
      PLIN; printf("%s: histogram %s not found\n",__FILE__,hNameH);
      return -3;
    }
    return 0;
}

void AliFastMuonTriggerEff::Init()
{
//
//  Initialization
//
    fZones=0;
    Char_t file[100]="$(ALICE_ROOT)/FASTSIM/data/MUONtriggerLUT_V2.4nvdn.root";
    printf("Initializing %s / %s\n", fName.Data(), fTitle.Data());
    printf("using data from file: %s\n",file);
    printf("AliFastMuonTriggerEff: Initialization with background level: %f\n",fBkg);
    TFile *f = new TFile(file);
    if(f->IsZombie()) {
        PLIN; printf("Cannot open file: %s\n",file);
        return;
    }
    f->ls();
    Int_t intb=0;
    Char_t namet[10];
    if(TMath::Abs(fBkg)<0.00001){
      sprintf(namet,"00");
    }else if(TMath::Abs(fBkg-0.5)<0.00001){
      sprintf(namet,"05");
    }else if(TMath::Abs(fBkg-1.0)<0.00001){
      sprintf(namet,"10");
    }else{
      PLIN; printf("A table for Bkg level: %f does not exists\n",fBkg);
      intb=1;
    }
    if(intb){ // Interpolation between background levels
      PLIN; printf("Interpolating Bkg level: %f\n",fBkg);
      TH3F* ha1,*hl1,*hh1,*ha2,*hl2,*hh2,*ha0,*hl0,*hh0;
      Char_t name1[10],name2[10]; Float_t b1,b2;
      if(fBkg>0&&fBkg<0.5){
        sprintf(name1,"00");
        sprintf(name2,"05");
	b1=0.;
	b2=0.5;
      }else if(fBkg>0.5){
        sprintf(name1,"05");
        sprintf(name2,"10");
	b1=0.5;
	b2=1.0;
	if(fBkg>1.0){
	  for(Int_t i=0; i<10;i++){
	    PLIN; printf("WARNING!!!! You are extrapolating above background 1.0\n");
	  }
	}
      }else{
        PLIN; printf("Bkg level: %f is not supported\n",fBkg);
        return;
      }
      if(LoadTables(name1)){
        PLIN; printf("Error in loading trigger tables\n");
	return;
      }
      PLIN; printf("We use tables for %f and %f to interpolate %f Bkg level\n",b1,b2,fBkg);
      ha0=(TH3F*)fhEffAPt->Clone("hEffAPtXX"); ha0->Reset();
      hl0=(TH3F*)fhEffLPt->Clone("hEffLPtXX"); hl0->Reset();
      hh0=(TH3F*)fhEffHPt->Clone("hEffHPtXX"); hh0->Reset();
      ha1=fhEffAPt;
      hl1=fhEffLPt;
      hh1=fhEffHPt;
      if(LoadTables(name2)){
        PLIN; printf("Error in loading trigger tables\n");
	return;
      }
      ha2=fhEffAPt;
      hl2=fhEffLPt;
      hh2=fhEffHPt;
      fhEffAPt=ha0;
      fhEffLPt=hl0;
      fhEffHPt=hh0;
      Int_t nnx=ha0->GetNbinsX()+1;
      Int_t nny=ha0->GetNbinsY()+1;
      Int_t nnz=ha0->GetNbinsZ()+1;
      for(Int_t ix=0; ix<=nnx; ix++){
	for(Int_t iy=0; iy<=nny; iy++){
	  for(Int_t iz=0; iz<=nnz; iz++){
	    Double_t y1,y2; Float_t cont;
	    y1=ha1->GetBinContent(ix,iy,iz); y2=ha2->GetBinContent(ix,iy,iz);
	    cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
	    fhEffAPt->SetBinContent(ix,iy,iz,cont);
	    y1=hl1->GetBinContent(ix,iy,iz); y2=hl2->GetBinContent(ix,iy,iz);
	    cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
	    fhEffLPt->SetBinContent(ix,iy,iz,cont);
	    y1=hh1->GetBinContent(ix,iy,iz); y2=hh2->GetBinContent(ix,iy,iz);
	    cont=Float_t(y1+(y2-y1)/(b2-b1)*(fBkg-b1)); if(cont>1)cont=1; if(cont<0)cont=0;
	    fhEffHPt->SetBinContent(ix,iy,iz,cont);
	  }
	}
      }
    }else{ // Use tables computed for selected backgound levels
      printf("Loading tables for background level: %f\n",fBkg);
      if(LoadTables(namet)){
        PLIN; printf("Error in loading trigger tables\n");
	return;
      }
    }
    fhEffAPt->SetDirectory(0);
    fhEffLPt->SetDirectory(0);
    fhEffHPt->SetDirectory(0);
    fhLX=fhEffLPt->GetXaxis();
    fhLY=fhEffLPt->GetYaxis();
    fhLZ=fhEffLPt->GetZaxis();
//
//
    if(f->Get("Description"))
    {
      fDescription=((TObjString*)f->Get("Description"))->GetString();
      printf("%s\n",fDescription.Data());
    }

    fThetaMin = fhEffLPt->GetXaxis()->GetXmin();
    fThetaMax = fhEffLPt->GetXaxis()->GetXmax();
    fnthetab=fhEffLPt->GetNbinsX();
    fDtheta   = (fThetaMax-fThetaMin)/fnthetab;

    fPhiMin   = fhEffLPt->GetYaxis()->GetXmin();
    fPhiMax   = fhEffLPt->GetYaxis()->GetXmax();
    fnphib=fhEffLPt->GetNbinsY();
    fDphi     = (fPhiMax-fPhiMin)/fnphib;

    fPtMin=fhEffLPt->GetZaxis()->GetXmin();
    fPtMax=fhEffLPt->GetZaxis()->GetXmax();
    fnptb=fhEffLPt->GetNbinsZ();
    fDpt      = (fPtMax-fPtMin)/fnptb;

    printf("%4d bins in theta [%f:%f]\n",fnthetab,fThetaMin,fThetaMax);
    printf("%4d bins in phi [%f:%f]\n",fnphib,fPhiMin,fPhiMax);
    printf("%4d bins in pt [%f:%f]\n",fnptb,fPtMin,fPtMax);

    fZones=fnthetab*fnphib;

    f->Close();
    if(fInt==0) {
      printf("Interpolation of trigger efficiencies is off!\n");
    } else {
      printf("Interpolation of trigger efficiencies is on!\n");
    }
}

void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,Float_t theta,
              Float_t phi, Float_t& effLow, Float_t& effHigh, Float_t& effAny)
{
    //
    //  Trigger efficiency for pt, theta, phi (low, high and "any" cut)
    //
#ifdef MYTRIGDEBUG
    printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
#endif
    effLow=0.;
    effHigh=0.;
    effAny=0;
    if(fZones==0) {
        printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
	return;
    }
    if(pt<0) {
        printf("Warning: pt: %f < 0. GeV/c\n",pt);
	return;	
    }

    Int_t iPt   = fhLZ->FindBin((Double_t)pt);
    if(iPt>fnptb)iPt=fnptb;
    Int_t iPhi  = Int_t((phi-fPhiMin)/fDphi);
    if(phi<fPhiMin)iPhi=iPhi-1;
    Int_t iTheta = fhLX->FindBin((Double_t)theta);
#ifdef MYTRIGDEBUG
    printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
    printf(" 0:%1d iPt iTheta iPhi: %d %d %d\n",fInt,iPt,iTheta,iPhi);
#endif
    iPhi=iPhi-2*fnphib*(iPhi/(2*fnphib));
#ifdef MYTRIGDEBUG
    printf(" 1:%1d iPhi converted to: %d for angle equivalence\n",fInt,iPhi);
#endif
    if(iPhi<0)iPhi=-iPhi-1;
    if(iPhi>(fnphib-1))iPhi=2*fnphib-1-iPhi;
#ifdef MYTRIGDEBUG
    printf(" 2:%1d iPhi converted to: %d for the symmetry of the spectrometer\n",fInt,iPhi);
#endif
    if(charge==1.){
    } else if(charge==-1.) {
    iPhi=fnphib-1-iPhi;
#ifdef MYTRIGDEBUG
    printf(" 3:%1d iPhi converted to: %d for the charge symmetry\n",fInt,iPhi);
#endif
    } else {
        printf("Warning: not understand charge: %f\n",charge);
        return;
    }
    if(iTheta<=0||iTheta>fnthetab) {
        printf("Warning: theta: %f outside acceptance\n",theta);
        return;
    }
    if(iPt<0) {
        printf("Warning: what do you mean with pt: %f <0?\n",pt);
	return;
    }
    iPhi++;
#ifdef MYTRIGDEBUG
    printf(" 4:%1d Getting: iTheta, iPhi, iPt: %d %d %d\n",
              fInt,iTheta,iPhi,iPt);
#endif
    effLow =fhEffLPt->GetBinContent(iTheta,iPhi,iPt);
    effHigh=fhEffHPt->GetBinContent(iTheta,iPhi,iPt);
    effAny =fhEffAPt->GetBinContent(iTheta,iPhi,iPt);
#ifdef MYTRIGDEBUG
    printf(" 4:%1d Result: charge, iTheta, iPhi, iPt: %f %d %d %d effLow: %f, effHigh: %f, effAny: %f\n",
              fInt,charge,iTheta,iPhi,iPt,effLow,effHigh,effAny);
#endif
        
    if(fInt==1) {
      Float_t angl,angh,anga;
      Float_t effLowp,effHighp,effAnyp;
      Float_t ptc=(iPt+0.5)*fDpt;  // The center of current bin
      #ifdef MYTRIGDEBUG
        printf(" 5:1 The center of current bin iPt: %d is: %f\n",iPt,ptc);
      #endif
      if(iPt==fnptb) {
        #ifdef MYTRIGDEBUG
	printf(" 6:1 No more points above! No interpolation is needed!\n");
        #endif
	return;        
      }else if(ptc==pt){
        #ifdef MYTRIGDEBUG
	  printf(" 6:1 No interpolation is needed!\n");
        #endif
	return;
      }else if(ptc>pt){
	// Looking for previous point
        if(iPt>1) {
          effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt-1);
          effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt-1);
          effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt-1);
          #ifdef MYTRIGDEBUG
	  printf(" 7:1 A simple look to previous point: %d: %f %f\n",iPt-1,effLowp,effHighp);
          #endif
	} else {
	  effLowp=0.;
	  effHighp=0.;
          effAnyp=0;
          #ifdef MYTRIGDEBUG
          printf(" 8:1 result is: %f %f %f\n",effLowp,effHighp,effAnyp);
	  #endif	  
	}
        angl=(effLow-effLowp)/fDpt;
        angh=(effHigh-effHighp)/fDpt;
        anga=(effAny-effAnyp)/fDpt;
      }else{
	// Looking for next point
	if(iPt<fnptb) {
          effLowp =fhEffLPt->GetBinContent(iTheta,iPhi,iPt+1);
          effHighp=fhEffHPt->GetBinContent(iTheta,iPhi,iPt+1);
          effAnyp =fhEffAPt->GetBinContent(iTheta,iPhi,iPt+1);
          #ifdef MYTRIGDEBUG
	  printf(" 7:1 A simple look to next point: %d: %f %f %f\n",iPt-1,effLowp,effHighp,effAnyp);
          #endif
	} else {
	  effLowp=effLow;
	  effHighp=effHigh;
	  effAnyp=effAny;
          #ifdef MYTRIGDEBUG
          printf(" 8:1 result is: pt: %f %f %f\n",effLowp,effHighp,effAnyp);
	  #endif	  
        }
        angl=(effLowp-effLow)/fDpt;
        angh=(effHighp-effHigh)/fDpt;
        anga=(effAnyp-effAny)/fDpt;
      }
      effLow=effLow+angl*(pt-ptc);
      effHigh=effHigh+angh*(pt-ptc);
      effAny=effAny+anga*(pt-ptc);
      #ifdef MYTRIGDEBUG
      printf(" 9:1 the interpolation coefficients are: %f %f %f\n",angl,angh,anga);
      #endif
  }
  #ifdef MYTRIGDEBUG
  printf("10:%1d effLow, effHigh=%f %f %f\n",fInt,effLow,effHigh,effAny);
  #endif
  return;
}


Float_t AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,
                   Float_t theta, Float_t phi)
{
    //
    // Trigger efficiency for pt, theta, phi depending of fCut
    // 
    if(fZones==0) {
        printf("Call to uninitialized object of class: AliFastMuonTriggerEff\n");
	return 0.;
    }
    Float_t eff;
    Float_t effLow, effHigh, effAny;
    
    Evaluate(charge,pt,theta,phi,effLow,effHigh,effAny);
    if (fCut == kLow) 
	eff  = effLow;
    else if (fCut == kHigh)
	eff  = effHigh;
    else if (fCut == kAny)
        eff  = effAny;
    else
        eff  = 0;

    return eff;
}

AliFastMuonTriggerEff& AliFastMuonTriggerEff::operator=(const  AliFastMuonTriggerEff& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return *this;
}
Commit	Line	Data
6255180c	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
a42548b0	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().
6255180c	26
490243f3	27
490243f3	28	#include "AliFastMuonTriggerEff.h"
2724a683	29	#include "TROOT.h"
	30	#include "TFile.h"
	31	#include "stdlib.h"
	32	#include "TH3.h"
	33	#include "TObjString.h"
	34
490243f3	35	#define PLIN printf("%s: %d: ",__FILE__,__LINE__)
6255180c	36
ecda8707	37	ClassImp(AliFastMuonTriggerEff)
6255180c	38
6255180c	39	AliFastMuonTriggerEff::AliFastMuonTriggerEff():
e6e76983	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)
6255180c	68	{
	69	//
	70	// Default constructor
	71	//
e6e76983	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
6255180c	105	}
6255180c	106
a42548b0	107	AliFastMuonTriggerEff::AliFastMuonTriggerEff(const AliFastMuonTriggerEff& eff)
e6e76983	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)
a42548b0	136	{
	137	// Copy constructor
	138	eff.Copy(*this);
	139	}
	140
490243f3	141	void AliFastMuonTriggerEff::SetCut(Int_t cut)
490243f3	142	{
2724a683	143	//
2724a683	144	// Set the pt cut
490243f3	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	}
6255180c	158
490243f3	159	Int_t AliFastMuonTriggerEff::SetBkgLevel(Float_t Bkg)
6255180c	160	{
2724a683	161	//
	162	// Set the background level
	163	//
490243f3	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	}
6255180c	181
490243f3	182	Int_t AliFastMuonTriggerEff::ForceBkgLevel(Float_t Bkg)
490243f3	183	{
2724a683	184	//
	185	// Check and enforce consistency of the background level
	186	//
490243f3	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	}
6255180c	202
490243f3	203	Int_t AliFastMuonTriggerEff::LoadTables(Char_t *namet=""){
2724a683	204	//
	205	// Load the trigger tables
	206	//
490243f3	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;
6255180c	227	}
6255180c	228
490243f3	229	void AliFastMuonTriggerEff::Init()
6255180c	230	{
6255180c	231	//
490243f3	232	// Initialization
6255180c	233	//
490243f3	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);
6255180c	239	TFile *f = new TFile(file);
6255180c	240	if(f->IsZombie()) {
490243f3	241	PLIN; printf("Cannot open file: %s\n",file);
6255180c	242	return;
	243	}
	244	f->ls();
490243f3	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;
6255180c	256	}
490243f3	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	}
6255180c	275	}
490243f3	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	}
6255180c	318	}
490243f3	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	}
6255180c	326	}
490243f3	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());
6255180c	339	}
490243f3	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
6255180c	362	f->Close();
490243f3	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	}
6255180c	368	}
6255180c	369
490243f3	370	void AliFastMuonTriggerEff::Evaluate(Float_t charge, Float_t pt,Float_t theta,
490243f3	371	Float_t phi, Float_t& effLow, Float_t& effHigh, Float_t& effAny)
6255180c	372	{
490243f3	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
6255180c	379	effLow=0.;
6255180c	380	effHigh=0.;
490243f3	381	effAny=0;
6255180c	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	}
490243f3	390
	391	Int_t iPt = fhLZ->FindBin((Double_t)pt);
	392	if(iPt>fnptb)iPt=fnptb;
6255180c	393	Int_t iPhi = Int_t((phi-fPhiMin)/fDphi);
6255180c	394	if(phi<fPhiMin)iPhi=iPhi-1;
490243f3	395	Int_t iTheta = fhLX->FindBin((Double_t)theta);
6255180c	396	#ifdef MYTRIGDEBUG
490243f3	397	printf("Evaluate(ch=%2.0f, pt=%10.6f, theta=%7.2f, phi=%8.2f ...)\n",charge,pt,theta,phi);
490243f3	398	printf(" 0:%1d iPt iTheta iPhi: %d %d %d\n",fInt,iPt,iTheta,iPhi);
6255180c	399	#endif
490243f3	400	iPhi=iPhi-2fnphib(iPhi/(2*fnphib));
6255180c	401	#ifdef MYTRIGDEBUG
490243f3	402	printf(" 1:%1d iPhi converted to: %d for angle equivalence\n",fInt,iPhi);
6255180c	403	#endif
6255180c	404	if(iPhi<0)iPhi=-iPhi-1;
490243f3	405	if(iPhi>(fnphib-1))iPhi=2*fnphib-1-iPhi;
6255180c	406	#ifdef MYTRIGDEBUG
490243f3	407	printf(" 2:%1d iPhi converted to: %d for the symmetry of the spectrometer\n",fInt,iPhi);
6255180c	408	#endif
	409	if(charge==1.){
	410	} else if(charge==-1.) {
490243f3	411	iPhi=fnphib-1-iPhi;
6255180c	412	#ifdef MYTRIGDEBUG
490243f3	413	printf(" 3:%1d iPhi converted to: %d for the charge symmetry\n",fInt,iPhi);
6255180c	414	#endif
	415	} else {
	416	printf("Warning: not understand charge: %f\n",charge);
	417	return;
	418	}
490243f3	419	if(iTheta<=0\|\|iTheta>fnthetab) {
6255180c	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	}
490243f3	427	iPhi++;
6255180c	428	#ifdef MYTRIGDEBUG
490243f3	429	printf(" 4:%1d Getting: iTheta, iPhi, iPt: %d %d %d\n",
490243f3	430	fInt,iTheta,iPhi,iPt);
6255180c	431	#endif
490243f3	432	effLow =fhEffLPt->GetBinContent(iTheta,iPhi,iPt);
	433	effHigh=fhEffHPt->GetBinContent(iTheta,iPhi,iPt);
	434	effAny =fhEffAPt->GetBinContent(iTheta,iPhi,iPt);
6255180c	435	#ifdef MYTRIGDEBUG
490243f3	436	printf(" 4:%1d Result: charge, iTheta, iPhi, iPt: %f %d %d %d effLow: %f, effHigh: %f, effAny: %f\n",
490243f3	437	fInt,charge,iTheta,iPhi,iPt,effLow,effHigh,effAny);
6255180c	438	#endif
490243f3	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
bc77d57e	511
6255180c	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;
490243f3	524	Float_t effLow, effHigh, effAny;
6255180c	525
490243f3	526	Evaluate(charge,pt,theta,phi,effLow,effHigh,effAny);
6255180c	527	if (fCut == kLow)
6255180c	528	eff = effLow;
490243f3	529	else if (fCut == kHigh)
6255180c	530	eff = effHigh;
490243f3	531	else if (fCut == kAny)
	532	eff = effAny;
	533	else
	534	eff = 0;
6255180c	535
	536	return eff;
	537	}
a42548b0	538
	539	AliFastMuonTriggerEff& AliFastMuonTriggerEff::operator=(const AliFastMuonTriggerEff& rhs)
	540	{
	541	// Assignment operator
	542	rhs.Copy(*this);
	543	return *this;
	544	}
	545