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