PWG2/FORWARD/analysis2/AliFMDCorrector.cxx

   1 //
   2 // This class calculates the exclusive charged particle density
   3 // in each for the 5 FMD rings.
   4 //
   5 #include "AliFMDCorrector.h"
   6 #include <AliESDFMD.h>
   7 #include <TAxis.h>
   8 #include <TList.h>
   9 #include <TMath.h>
  10 #include "AliForwardCorrectionManager.h"
  11 #include "AliLog.h"
  12 #include <TH2D.h>
  13 #include <TROOT.h>
  14 #include <iostream>
  15 #include <iomanip>
  16
  17 ClassImp(AliFMDCorrector)
  18 #if 0
  19 ; // For Emacs
  20 #endif
  21
  22 //____________________________________________________________________
  23 AliFMDCorrector::AliFMDCorrector()
  24   : TNamed(),
  25     fRingHistos(),
  26     fUseMergingEfficiency(true),
  27     fDebug(0)
  28 {
  29   // Constructor
  30 }
  31
  32 //____________________________________________________________________
  33 AliFMDCorrector::AliFMDCorrector(const char* title)
  34   : TNamed("fmdCorrector", title),
  35     fRingHistos(),
  36     fUseMergingEfficiency(true),
  37     fDebug(0)
  38 {
  39   // Constructor
  40   //
  41   // Parameters:
  42   //   title   Title
  43   fRingHistos.SetName(GetName());
  44   fRingHistos.Add(new RingHistos(1, 'I'));
  45   fRingHistos.Add(new RingHistos(2, 'I'));
  46   fRingHistos.Add(new RingHistos(2, 'O'));
  47   fRingHistos.Add(new RingHistos(3, 'I'));
  48   fRingHistos.Add(new RingHistos(3, 'O'));
  49 }
  50
  51 //____________________________________________________________________
  52 AliFMDCorrector::AliFMDCorrector(const AliFMDCorrector& o)
  53   : TNamed(o),
  54     fRingHistos(),
  55     fUseMergingEfficiency(o.fUseMergingEfficiency),
  56     fDebug(o.fDebug)
  57 {
  58   // Copy constructor
  59   //
  60   // Parameters:
  61   //  o  Object to copy from
  62   TIter    next(&o.fRingHistos);
  63   TObject* obj = 0;
  64   while ((obj = next())) fRingHistos.Add(obj);
  65 }
  66
  67 //____________________________________________________________________
  68 AliFMDCorrector::~AliFMDCorrector()
  69 {
  70   // Destructor
  71   //
  72   //
  73   fRingHistos.Delete();
  74 }
  75
  76 //____________________________________________________________________
  77 AliFMDCorrector&
  78 AliFMDCorrector::operator=(const AliFMDCorrector& o)
  79 {
  80   // Assignment operator
  81   //
  82   // Parameters:
  83   //   o   Object to assign from
  84   TNamed::operator=(o);
  85
  86   fDebug   = o.fDebug;
  87   fRingHistos.Delete();
  88   fUseMergingEfficiency = o.fUseMergingEfficiency;
  89   TIter    next(&o.fRingHistos);
  90   TObject* obj = 0;
  91   while ((obj = next())) fRingHistos.Add(obj);
  92
  93   return *this;
  94 }
  95
  96 //____________________________________________________________________
  97 AliFMDCorrector::RingHistos*
  98 AliFMDCorrector::GetRingHistos(UShort_t d, Char_t r) const
  99 {
 100   //
 101   // Get the ring histogram container
 102   // Parameters:
 103   //    d Detector
 104   //    r Ring
 105   //
 106   // Return:
 107   //    Ring histogram container
 108   //
 109   Int_t idx = -1;
 110   switch (d) {
 111   case 1: idx = 0; break;
 112   case 2: idx = 1 + (r == 'I' || r == 'i' ? 0 : 1); break;
 113   case 3: idx = 3 + (r == 'I' || r == 'i' ? 0 : 1); break;
 114   }
 115   if (idx < 0 || idx >= fRingHistos.GetEntries()) return 0;
 116
 117   return static_cast<RingHistos*>(fRingHistos.At(idx));
 118 }
 119
 120 //____________________________________________________________________
 121 Bool_t
 122 AliFMDCorrector::Correct(AliForwardUtil::Histos& hists,
 123                            UShort_t                vtxbin)
 124 {
 125   //
 126   // Do the calculations
 127   // Parameters:
 128   //    hists    Cache of histograms
 129   //    vtxBin   Vertex bin
 130   //
 131   // Return:
 132   //    true on successs
 133   //
 134   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
 135
 136   UShort_t uvb = vtxbin;
 137   for (UShort_t d=1; d<=3; d++) {
 138     UShort_t nr = (d == 1 ? 1 : 2);
 139     for (UShort_t q=0; q<nr; q++) {
 140       Char_t      r  = (q == 0 ? 'I' : 'O');
 141       TH2D*       h  = hists.Get(d,r);
 142       RingHistos* rh = GetRingHistos(d,r);
 143       TH2D*       bg = fcm.GetSecondaryMap()->GetCorrection(d, r, uvb);
 144       TH2D*       ef = fcm.GetVertexBias()->GetCorrection(r, uvb);
 145       TH2D*       ac = fcm.GetAcceptance()->GetCorrection(d, r, uvb);
 146       if (!bg) {
 147         AliWarning(Form("No secondary correction for FMDM%d%c in vertex bin %d",
 148                         d, r, uvb));
 149         continue;
 150       }
 151       if (!ef) {
 152         AliWarning(Form("No event %s vertex bias correction in vertex bin %d",
 153                         (r == 'I' || r == 'i' ? "inner" : "outer"), uvb));
 154         continue;
 155       }
 156       if (!ac) {
 157         AliWarning(Form("No acceptance correction for FMD%d%c in vertex bin %d",
 158                         d, r, uvb));
 159         continue;
 160       }
 161
 162       // Divide by primary/total ratio
 163       h->Divide(bg);
 164
 165       // Divide by the event selection efficiency
 166       h->Divide(ef);
 167
 168       // Divide by the acceptance correction
 169       h->Divide(ac);
 170
 171       if (fUseMergingEfficiency) {
 172         if (!fcm.GetMergingEfficiency()) {
 173           AliWarning("No merging efficiencies");
 174           continue;
 175         }
 176         TH1D* sf = fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb);
 177         if (!fcm.GetMergingEfficiency()->GetCorrection(d,r,uvb)) {
 178           AliWarning(Form("No merging efficiency for FMD%d%c at vertex bin %d",
 179                           d, r, uvb));
 180           continue;
 181         }
 182
 183
 184         for (Int_t ieta = 1; ieta <= h->GetNbinsX(); ieta++) {
 185           Float_t c  = sf->GetBinContent(ieta);
 186           Float_t ec = sf->GetBinError(ieta);
 187
 188           if (c == 0) continue;
 189
 190           for (Int_t iphi = 1; iphi <= h->GetNbinsY(); iphi++) {
 191             Double_t m  = h->GetBinContent(ieta, iphi) / c;
 192             Double_t em = h->GetBinError(ieta, iphi);
 193
 194             Double_t e  = TMath::Sqrt(em * em + (m * ec) * (m * ec));
 195
 196             h->SetBinContent(ieta,iphi,m);
 197             h->SetBinError(ieta,iphi,e);
 198           }
 199         }
 200       }
 201       rh->fDensity->Add(h);
 202     }
 203   }
 204
 205   return kTRUE;
 206 }
 207
 208 //____________________________________________________________________
 209 void
 210 AliFMDCorrector::ScaleHistograms(TList* dir, Int_t nEvents)
 211 {
 212   //
 213   // Scale the histograms to the total number of events
 214   // Parameters:
 215   //    dir     Where the output is stored
 216   //    nEvents Number of events
 217   //
 218   if (nEvents <= 0) return;
 219   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 220   if (!d) return;
 221
 222   TIter    next(&fRingHistos);
 223   RingHistos* o = 0;
 224   while ((o = static_cast<RingHistos*>(next())))
 225     o->ScaleHistograms(d, nEvents);
 226 }
 227 //____________________________________________________________________
 228 void
 229 AliFMDCorrector::DefineOutput(TList* dir)
 230 {
 231   TList* d = new TList;
 232   d->SetName(GetName());
 233   dir->Add(d);
 234   TIter    next(&fRingHistos);
 235   RingHistos* o = 0;
 236   while ((o = static_cast<RingHistos*>(next()))) {
 237     o->Output(d);
 238   }
 239 }
 240
 241 //____________________________________________________________________
 242 void
 243 AliFMDCorrector::Print(Option_t* /* option */) const
 244 {
 245   //
 246   // Print information
 247   // Parameters:
 248   //    option Not used
 249   //
 250   char ind[gROOT->GetDirLevel()+1];
 251   for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
 252   ind[gROOT->GetDirLevel()] = '\0';
 253   std::cout << ind << "AliFMDCorrector: " << GetName() <<  std::endl;
 254 }
 255
 256 //====================================================================
 257 AliFMDCorrector::RingHistos::RingHistos()
 258   : AliForwardUtil::RingHistos(),
 259     fDensity(0)
 260 {
 261   // Constructor
 262   //
 263   //
 264 }
 265
 266 //____________________________________________________________________
 267 AliFMDCorrector::RingHistos::RingHistos(UShort_t d, Char_t r)
 268   : AliForwardUtil::RingHistos(d,r),
 269     fDensity(0)
 270 {
 271   //
 272   // Constructor
 273   // Parameters:
 274   //    d detector
 275   //    r ring
 276   //
 277   fDensity = new TH2D(Form("FMD%d%c_Primary_Density", d, r),
 278                       Form("in FMD%d%c", d, r),
 279                       200, -4, 6, (r == 'I' || r == 'i' ? 20 : 40),
 280                       0, 2*TMath::Pi());
 281   fDensity->SetDirectory(0);
 282   fDensity->SetXTitle("#eta");
 283   fDensity->SetYTitle("#phi [radians]");
 284   fDensity->SetZTitle("Primary N_{ch} density");
 285 }
 286 //____________________________________________________________________
 287 AliFMDCorrector::RingHistos::RingHistos(const RingHistos& o)
 288   : AliForwardUtil::RingHistos(o),
 289     fDensity(o.fDensity)
 290 {
 291   //
 292   // Copy constructor
 293   // Parameters:
 294   //    o Object to copy from
 295   //
 296 }
 297
 298 //____________________________________________________________________
 299 AliFMDCorrector::RingHistos&
 300 AliFMDCorrector::RingHistos::operator=(const RingHistos& o)
 301 {
 302   //
 303   // Assignment operator
 304   // Parameters:
 305   //    o Object to assign from
 306   //
 307   // Return:
 308   //    Reference to this
 309   //
 310   AliForwardUtil::RingHistos::operator=(o);
 311
 312   if (fDensity) delete fDensity;
 313
 314   fDensity = static_cast<TH2D*>(o.fDensity->Clone());
 315
 316   return *this;
 317 }
 318 //____________________________________________________________________
 319 AliFMDCorrector::RingHistos::~RingHistos()
 320 {
 321   //
 322   // Destructor
 323   //
 324   if (fDensity) delete fDensity;
 325 }
 326
 327 //____________________________________________________________________
 328 void
 329 AliFMDCorrector::RingHistos::Output(TList* dir)
 330 {
 331   //
 332   // Make output
 333   // Parameters:
 334   //    dir Where to put it
 335   //
 336   TList* d = DefineOutputList(dir);
 337   d->Add(fDensity);
 338 }
 339
 340 //____________________________________________________________________
 341 void
 342 AliFMDCorrector::RingHistos::ScaleHistograms(TList* dir, Int_t nEvents)
 343 {
 344   //
 345   // Scale the histograms to the total number of events
 346   // Parameters:
 347   //    dir     where the output is stored
 348   //    nEvents Number of events
 349   //
 350   TList* l = GetOutputList(dir);
 351   if (!l) return;
 352
 353   TH1* density = GetOutputHist(l,Form("%s_Primary_Density", fName.Data()));
 354   if (density) density->Scale(1./nEvents);
 355 }
 356
 357 //____________________________________________________________________
 358 //
 359 // EOF
 360 //
 361
 362
 363