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