PWGDQ/base/AliQuarkoniaEfficiency.cxx

   1
   2 /* $Id$ */
   3
   4 //===================================================================
   5 //  Class AliQuarkoniaEfficiency
   6 //
   7 //  This class will provide the quarkonia reconstruction efficiency
   8 //  in ALICE without acceptance consideration.
   9 //
  10 //
  11 //  Reconstruction efficiency has been evaluated by means of a flat
  12 //  y and pt distribution of quarkonia in -4 < y < -2.5,
  13 //  0 < pt < 20 GeV/c. Weights have been used to evaluate the
  14 //  reconstruction efficiency in different parameterizations.
  15 //
  16 //  Beware that efficiency histos are preliminary.
  17 //  Just Jpsi, Dimuon, UnlikePair trigger
  18 //  efficiencies should be considered.
  19 //
  20 //
  21 //  Example:
  22 //   Double_t eff,error;
  23 //   AliQuarkoniaEfficiency * JPsiEff = new AliQuarkoniaEfficiency();
  24 //   JPsiEff->Init();
  25 //   JPsiEff->GetEfficiency(Rapidity,Pt,eff,error);
  26 //   printf(" eff = %2.2e  error %2.2e \n",eff,error);
  27 //
  28 //
  29 //  Subatech 2006
  30 //===================================================================
  31
  32 // Root
  33 #include "TFile.h"
  34 #include "TH2.h"
  35
  36 // AliRoot includes
  37 #include "AliQuarkoniaEfficiency.h"
  38 #include "AliLog.h"
  39
  40
  41 //_______________________________________________________________________
  42 AliQuarkoniaEfficiency::AliQuarkoniaEfficiency(Int_t quarkoniaResonance, Int_t decayChannel,
  43                                                Int_t simParameterization):
  44   fEfficiencyFileName("QuarkoniaEfficiency.root"),
  45   fQuarkoniaResonance(quarkoniaResonance),
  46   fDecayChannel(decayChannel),
  47   fParameterization(simParameterization),
  48   fTriggerType(kPairUnlikeApt),
  49   fTrigger(kFALSE),
  50   fEfficiency(0x0)
  51 {
  52   // Constructor
  53 }
  54
  55 //_______________________________________________________________________
  56 AliQuarkoniaEfficiency::~AliQuarkoniaEfficiency()
  57 {
  58   // Destructor
  59   delete fEfficiency;
  60 }
  61
  62 //_______________________________________________________________________
  63 void AliQuarkoniaEfficiency::Init()
  64 {
  65   // Initialize method
  66   switch (fQuarkoniaResonance)  {
  67   case kJpsi:
  68     SetTitle("Jpsi");
  69     break;
  70   case kPsiP:
  71     SetTitle("PsiP");
  72     break;
  73   case kUpsilon:
  74     SetTitle("Upsilon");
  75     break;
  76   case kUpsilonP:
  77     SetTitle("UpsilonP");
  78     break;
  79   case kUpsilonPP:
  80     SetTitle("UpsilonPP");
  81     break;
  82   case kOmega:
  83     SetTitle("Omega");
  84     break;
  85   case kPhi:
  86     SetTitle("Phi");
  87     break;
  88   }
  89
  90   switch (fDecayChannel) {
  91   case kDimuon:
  92     SetName("Dimuon");
  93     break;
  94   case kDielectron:
  95     SetName("Dielectron");
  96     break;
  97   }
  98
  99   const char *param=0;
 100   switch (fParameterization){
 101   case kFlat:
 102     param = "Flat";
 103     break;
 104   case kCDFscaled:
 105     param = "CDFscaled";
 106     break;
 107   case kCDFscaledPP:
 108     param = "CDFscaledPP";
 109     break;
 110   }
 111
 112   const char *trig=0;
 113   switch (fTriggerType){
 114   case kSinglePlusLpt:
 115     trig = "SinglePlusLpt";
 116     break;
 117   case kSinglePlusHpt:
 118     trig = "SinglePlusHpt";
 119     break;
 120   case kSinglePlusApt:
 121     trig = "SinglePlusApt";
 122     break;
 123   case kSingleMinusLpt:
 124     trig = "SingleMinusLpt";
 125     break;
 126   case kSingleMinusHpt:
 127     trig = "SingleMinusHpt";
 128     break;
 129   case kSingleMinusApt:
 130     trig = "SingleMinusApt";
 131     break;
 132   case kSingleUndefLpt:
 133     trig = "SingleUndefLpt";
 134     break;
 135   case kSingleUndefHpt:
 136     trig = "SingleUndefHpt";
 137     break;
 138   case kSingleUndefApt:
 139     trig = "SingleUndefApt";
 140     break;
 141   case kPairUnlikeLpt:
 142     trig = "PairUnlikeLpt";
 143     break;
 144   case kPairUnlikeHpt:
 145     trig = "PairUnlikeHpt";
 146     break;
 147   case kPairUnlikeApt:
 148     trig = "PairUnlikeApt";
 149     break;
 150   case kPairLikeLpt:
 151     trig = "PairLikeLpt";
 152     break;
 153   case kPairLikeHpt:
 154     trig = "PairLikeHpt";
 155     break;
 156   case kPairLikeApt:
 157     trig = "PairLikeApt";
 158     break;
 159   }
 160
 161
 162   if(!fEfficiency) delete fEfficiency;
 163
 164   TFile efficiencyFile(fEfficiencyFileName);
 165   if ( efficiencyFile.IsOpen() ) {
 166
 167     char quarkoniaDir[15];
 168     snprintf(quarkoniaDir,15,"%s",GetTitle());
 169     if (! efficiencyFile.cd(quarkoniaDir) ){
 170       AliError(Form("Directory %s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found ",
 171                     quarkoniaDir,fEfficiencyFileName.Data(),GetTitle(),GetName() ));
 172       return;
 173     }
 174
 175     char histosDir[30];
 176     snprintf(histosDir,30,"%s/%s_%s_%s",quarkoniaDir,GetTitle(),GetName(),param);
 177     if(! efficiencyFile.cd(histosDir) ){
 178       AliError(Form("Subdirectory %s/%s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found ",
 179                     quarkoniaDir,histosDir,fEfficiencyFileName.Data(),GetTitle(),GetName() ));
 180       return;
 181     }
 182
 183     char histoname[50];
 184     if(fTrigger) snprintf(histoname,50,"h%sEfficiencyPtRap_%s",GetTitle(),trig);
 185     else snprintf(histoname,50,"h%sEfficiencyPtRap",GetTitle());
 186     char histonameposition[99];
 187     snprintf(histonameposition,99,"%s/%s",histosDir,histoname);
 188     fEfficiency = (TH2F*)efficiencyFile.Get(histonameposition);
 189
 190     if ( !fEfficiency ) {
 191       AliError(Form("Histo %s not found in file %s \n Efficiency data for quarkonia %s and channel %s not found",
 192                     histoname, fEfficiencyFileName.Data(), GetTitle(), GetName() ));
 193     }
 194     else {
 195       fEfficiency->SetDirectory(0);
 196     }
 197     efficiencyFile.Close();
 198
 199   }
 200   else {
 201     AliError(Form("File %s not found",fEfficiencyFileName.Data()));
 202   }
 203
 204 }
 205
 206 //_______________________________________________________________________
 207 TH2F*  AliQuarkoniaEfficiency::GetEfficiencyHisto() const
 208 {
 209   // Returns the efficiency histogram
 210   if (fEfficiency) return fEfficiency;
 211   else {
 212     AliError(Form("Efficiency data for quarkonia %s and channel %s not found",GetTitle(),GetName()));
 213     return 0x0;
 214   }
 215 }
 216
 217 //_______________________________________________________________________
 218 void  AliQuarkoniaEfficiency::GetEfficiency(Float_t rap, Float_t pT, Double_t &eff, Double_t &error)
 219 {
 220   // Evaluates the efficiency for a given (y,Pt) of the quarkonia
 221   Int_t binx=0;
 222   Int_t biny=0;
 223
 224   if (!fEfficiency) {
 225     AliError(Form("Efficiency data for quarkonia %s and channel %s not found",GetTitle(),GetName()));
 226   }
 227   else {
 228     if ( rap < (fEfficiency->GetXaxis())->GetXmin()  ||
 229          rap > (fEfficiency->GetXaxis())->GetXmax()  ||
 230          pT  < (fEfficiency->GetYaxis())->GetXmin()  ||
 231          pT  > (fEfficiency->GetYaxis())->GetXmax()   ) {
 232       AliInfo("Values out of range");
 233       eff   = 0.;
 234       error = 0.;
 235     }
 236     else  {
 237       binx  = fEfficiency->GetXaxis()->FindBin(rap);
 238       biny  = fEfficiency->GetYaxis()->FindBin(pT);
 239       eff   = fEfficiency->GetBinContent(binx,biny);
 240       error = fEfficiency->GetBinError(binx,biny);
 241     }
 242   }
 243 }