PWGLF/FORWARD/analysis2/AliFMDMCCorrector.cxx

   1 //
   2 // This class calculates the exclusive charged particle density
   3 // in each for the 5 FMD rings.
   4 //
   5 // Input:
   6 //   - 5 RingHistos objects - each with a number of vertex dependent
   7 //     2D histograms of the inclusive charge particle density
   8 //
   9 // Output:
  10 //   - 5 RingHistos objects - each with a number of vertex dependent
  11 //     2D histograms of the exclusive charge particle density
  12 //
  13 // Corrections used:
  14 //   - AliFMDCorrSecondaryMap;
  15 //   - AliFMDCorrVertexBias
  16 //   - AliFMDCorrMergingEfficiency
  17 //
  18 #include "AliFMDMCCorrector.h"
  19 #include <AliESDFMD.h>
  20 #include <TAxis.h>
  21 #include <TList.h>
  22 #include <TMath.h>
  23 #include "AliForwardCorrectionManager.h"
  24 #include "AliFMDCorrSecondaryMap.h"
  25 #include "AliFMDCorrVertexBias.h"
  26 #include "AliLog.h"
  27 #include <TH2D.h>
  28 #include <TROOT.h>
  29 #include <TProfile2D.h>
  30 #include <iostream>
  31 ClassImp(AliFMDMCCorrector)
  32 #if 0
  33 ; // For Emacs
  34 #endif
  35
  36
  37 //____________________________________________________________________
  38 AliFMDMCCorrector::~AliFMDMCCorrector()
  39 {
  40   //
  41   // Destructor
  42   //
  43   // if (fComps) fComps->Clear();
  44   // if (fFMD1i) delete fFMD1i;
  45   // if (fFMD2i) delete fFMD2i;
  46   // if (fFMD2o) delete fFMD2o;
  47   // if (fFMD3i) delete fFMD3i;
  48   // if (fFMD3o) delete fFMD3o;
  49 }
  50
  51 //____________________________________________________________________
  52 AliFMDMCCorrector&
  53 AliFMDMCCorrector::operator=(const AliFMDMCCorrector& o)
  54 {
  55   //
  56   // Assignement operator
  57   //
  58   // Parameters:
  59   //    o Object to assign from
  60   //
  61   // Return:
  62   //    Reference to this object
  63   //
  64   if (&o == this) return *this;
  65   AliFMDCorrector::operator=(o);
  66   fSecondaryForMC = o.fSecondaryForMC;
  67   return *this;
  68 }
  69
  70 //____________________________________________________________________
  71 Bool_t
  72 AliFMDMCCorrector::CorrectMC(AliForwardUtil::Histos& hists,
  73                                UShort_t                vtxbin)
  74 {
  75   //
  76   // Do the calculations
  77   //
  78   // Parameters:
  79   //    hists    Cache of histograms
  80   //    vtxBin   Vertex bin
  81   //
  82   // Return:
  83   //    true on successs
  84   //
  85   if ((!fUseSecondaryMap || !fSecondaryForMC) && fUseVertexBias)
  86     return kTRUE;
  87
  88   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
  89
  90   UShort_t uvb = vtxbin;
  91   for (UShort_t d=1; d<=3; d++) {
  92     UShort_t nr = (d == 1 ? 1 : 2);
  93     for (UShort_t q=0; q<nr; q++) {
  94       Char_t      r  = (q == 0 ? 'I' : 'O');
  95       TH2D*       h  = hists.Get(d,r);
  96
  97
  98       if (fUseSecondaryMap && fSecondaryForMC) {
  99         TH2D*       bg = fcm.GetSecondaryMap()->GetCorrection(d,r,uvb);
 100         if (!bg) {
 101           AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
 102               d, r, uvb));
 103           continue;
 104         }
 105         // Divide by primary/total ratio
 106         h->Divide(bg);
 107       }
 108       if (fUseVertexBias) {
 109         TH2D*       ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
 110         if (!ef) {
 111           AliWarning(Form("No event vertex bias correction in vertex bin %d",
 112                         uvb));
 113           continue;
 114         }
 115         // Divide by the event selection efficiency
 116         h->Divide(ef);
 117       }
 118     }
 119   }
 120
 121   return kTRUE;
 122 }
 123
 124 //____________________________________________________________________
 125 void
 126 AliFMDMCCorrector::SetupForData(const TAxis& eAxis)
 127 {
 128   //
 129   // Initialize this object
 130   //
 131   // Parameters:
 132   //    etaAxis Eta axis to use
 133   //
 134   AliFMDCorrector::SetupForData(eAxis);
 135
 136   fFMD1i = Make(1,'I',eAxis);
 137   fFMD2i = Make(2,'I',eAxis);
 138   fFMD2o = Make(2,'O',eAxis);
 139   fFMD3i = Make(3,'I',eAxis);
 140   fFMD3o = Make(3,'O',eAxis);
 141
 142   fComps->Add(fFMD1i);
 143   fComps->Add(fFMD2i);
 144   fComps->Add(fFMD2o);
 145   fComps->Add(fFMD3i);
 146   fComps->Add(fFMD3o);
 147 }
 148
 149 //____________________________________________________________________
 150 TProfile2D*
 151 AliFMDMCCorrector::Make(UShort_t d, Char_t r,
 152                                 const TAxis& axis) const
 153 {
 154   //
 155   // MAke comparison profiles
 156   //
 157   // Parameters:
 158   //    d     Detector
 159   //    r     Ring
 160   //    axis  Eta axis
 161   //
 162   // Return:
 163   //    Newly allocated profile object
 164   //
 165   TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
 166                                    Form("ESD/MC signal for FMD%d%c", d, r),
 167                                    axis.GetNbins(),
 168                                    axis.GetXmin(),
 169                                    axis.GetXmax(),
 170                                    (r == 'I' || r == 'i') ? 20 : 40,
 171                                    0, 2*TMath::Pi());
 172   ret->GetXaxis()->SetTitle("#eta");
 173   ret->GetYaxis()->SetTitle("#varphi [degrees]");
 174   ret->GetZaxis()->SetTitle("#LT primary density ESD/MC#GT");
 175   ret->SetDirectory(0);
 176   return ret;
 177 }
 178 //____________________________________________________________________
 179 void
 180 AliFMDMCCorrector::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
 181 {
 182   //
 183   // Fill comparison profiles
 184   //
 185   // Parameters:
 186   //    d    Detector
 187   //    r    Ring
 188   //    esd  ESD histogram
 189   //    mc   MC histogram
 190   //
 191   if (!esd || !mc) return;
 192   TProfile2D* p = 0;
 193   switch (d) {
 194   case 1:  p = fFMD1i;                                   break;
 195   case 2:  p = (r == 'I' || r == 'i' ? fFMD2i : fFMD2o); break;
 196   case 3:  p = (r == 'I' || r == 'i' ? fFMD3i : fFMD3o); break;
 197   }
 198   if (!p) return;
 199
 200   for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
 201     Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
 202     for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
 203       Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
 204       Double_t mEsd = esd->GetBinContent(iEta,iPhi);
 205       Double_t mMc  = mc->GetBinContent(iEta,iPhi);
 206
 207       p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
 208     }
 209   }
 210 }
 211
 212 //____________________________________________________________________
 213 Bool_t
 214 AliFMDMCCorrector::CompareResults(AliForwardUtil::Histos& esd,
 215                                           AliForwardUtil::Histos& mc)
 216 {
 217   //
 218   // Compare the result of analysing the ESD for
 219   // the inclusive charged particle density to analysing
 220   // MC truth
 221   //
 222   // Parameters:
 223   //    esd
 224   //    mc
 225   //
 226   // Return:
 227   //   true
 228   //
 229   Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
 230   Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
 231   Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
 232   Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
 233   Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
 234
 235   return kTRUE;
 236 }
 237
 238 //____________________________________________________________________
 239 void
 240 AliFMDMCCorrector::CreateOutputObjects(TList* dir)
 241 {
 242   //
 243   // Output diagnostic histograms to directory
 244   //
 245   // Parameters:
 246   //    dir List to write in
 247   //
 248   AliFMDCorrector::CreateOutputObjects(dir);
 249   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 250
 251   fComps = new TList;
 252   fComps->SetName("esd_mc_comparison");
 253   d->Add(fComps);
 254 }
 255 #define PFB(N,FLAG)                             \
 256   do {                                                                  \
 257     AliForwardUtil::PrintName(N);                                       \
 258     std::cout << std::boolalpha << (FLAG) << std::noboolalpha << std::endl; \
 259   } while(false)
 260
 261 //____________________________________________________________________
 262 void
 263 AliFMDMCCorrector::Print(Option_t* option) const
 264 {
 265   //
 266   // Print information
 267   // Parameters:
 268   //    option Not used
 269   //
 270   AliFMDCorrector::Print(option);
 271   gROOT->IncreaseDirLevel();
 272   PFB("Use sec. map on MC", fSecondaryForMC);
 273   gROOT->DecreaseDirLevel();
 274 }
 275
 276 //____________________________________________________________________
 277 //
 278 // EOF
 279 //
 280
 281
 282