PWGLF/FORWARD/analysis2/AliFMDMCDensityCalculator.cxx

   1 //
   2 // This class calculates the inclusive charged particle density
   3 // in each for the 5 FMD rings based on the MC truth.
   4 //
   5 // Input:
   6 //   - AliMCEvent  MC truth event infromation
   7 //
   8 // Output:
   9 //   - None
  10 //
  11 // Corrections used:
  12 //   - None
  13 //
  14 //
  15 #include "AliFMDMCDensityCalculator.h"
  16 #include <TMath.h>
  17 #include "AliForwardCorrectionManager.h"
  18 #include "AliFMDStripIndex.h"
  19 #include "AliMCEvent.h"
  20 #include "AliESDFMD.h"
  21 #include "AliLog.h"
  22 #include <TH2D.h>
  23 #include <TProfile2D.h>
  24
  25 ClassImp(AliFMDMCDensityCalculator)
  26 #if 0
  27 ; // For Emacs
  28 #endif
  29
  30 //____________________________________________________________________
  31 AliFMDMCDensityCalculator::~AliFMDMCDensityCalculator()
  32 {
  33   //
  34   // Destructor
  35   //
  36   // if (fComps)  fComps->Delete();
  37 }
  38
  39 //____________________________________________________________________
  40 AliFMDMCDensityCalculator&
  41 AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
  42 {
  43   //
  44   // Assignement operator
  45   //
  46   // Parameters:
  47   //    o Object to assign from
  48   //
  49   // Return:
  50   //    Reference to this object
  51   //
  52   AliFMDDensityCalculator::operator=(o);
  53   return *this;
  54 }
  55
  56
  57 //____________________________________________________________________
  58 void
  59 AliFMDMCDensityCalculator::SetupForData(const TAxis& eAxis)
  60 {
  61   //
  62   // Initialize this object
  63   //
  64   // Parameters:
  65   //    etaAxis Eta axis to use
  66   //
  67   AliFMDDensityCalculator::SetupForData(eAxis);
  68   fFMD1i  = Make(1,'I',eAxis);
  69   fFMD2i  = Make(2,'I',eAxis);
  70   fFMD2o  = Make(2,'O',eAxis);
  71   fFMD3i  = Make(3,'I',eAxis);
  72   fFMD3o  = Make(3,'O',eAxis);
  73   fFMD1iC = Make(1,'I');
  74   fFMD2iC = Make(2,'I');
  75   fFMD2oC = Make(2,'O');
  76   fFMD3iC = Make(3,'I');
  77   fFMD3oC = Make(3,'O');
  78   fFMD1iD = Make(1,'I', 20);
  79   fFMD2iD = Make(2,'I', 20);
  80   fFMD2oD = Make(2,'O', 20);
  81   fFMD3iD = Make(3,'I', 20);
  82   fFMD3oD = Make(3,'O', 20);
  83
  84   fComps->Add(fFMD1i);
  85   fComps->Add(fFMD2i);
  86   fComps->Add(fFMD2o);
  87   fComps->Add(fFMD3i);
  88   fComps->Add(fFMD3o);
  89   fComps->Add(fFMD1iC);
  90   fComps->Add(fFMD2iC);
  91   fComps->Add(fFMD2oC);
  92   fComps->Add(fFMD3iC);
  93   fComps->Add(fFMD3oC);
  94   fComps->Add(fFMD1iD);
  95   fComps->Add(fFMD2iD);
  96   fComps->Add(fFMD2oD);
  97   fComps->Add(fFMD3iD);
  98   fComps->Add(fFMD3oD);
  99 }
 100
 101
 102 //____________________________________________________________________
 103 Bool_t
 104 AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
 105                                        AliForwardUtil::Histos& hists)
 106 {
 107   //
 108   // Calculate the charged particle density from the MC track references.
 109   //
 110   // Parameters:
 111   //    fmd    Forward MC event
 112   //    hists  Histograms to fill
 113   //    vz     Interaction z coordinate @f$ v_z@f$
 114   //    vtxBin bin corresponding to @f$ v_z@f$
 115   //
 116   // Return:
 117   //    true on success
 118   //
 119   for (UShort_t d=1; d<=3; d++) {
 120     UShort_t nr = (d == 1 ? 1 : 2);
 121     for (UShort_t q=0; q<nr; q++) {
 122       Char_t      r = (q == 0 ? 'I' : 'O');
 123       UShort_t    ns= (q == 0 ?  20 :  40);
 124       UShort_t    nt= (q == 0 ? 512 : 256);
 125       TH2D*       h = hists.Get(d,r);
 126
 127       for (UShort_t s=0; s<ns; s++) {
 128         for (UShort_t t=0; t<nt; t++) {
 129           Float_t mult = fmd.Multiplicity(d,r,s,t);
 130
 131           if (mult == 0 || mult > 20) continue;
 132
 133           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
 134           Float_t eta = fmd.Eta(d,r,s,t);
 135
 136           Float_t corr = Correction(d, r, t, eta, false);
 137           Float_t sig  = (corr <= 0 ? 0 : mult / corr);
 138           h->Fill(eta,phi,sig);
 139         }
 140       }
 141     }
 142   }
 143   return kTRUE;
 144 }
 145
 146 //____________________________________________________________________
 147 TProfile2D*
 148 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r,
 149                                 const TAxis& axis) const
 150 {
 151   //
 152   // MAke comparison profiles
 153   //
 154   // Parameters:
 155   //    d     Detector
 156   //    r     Ring
 157   //    axis  Eta axis
 158   //
 159   // Return:
 160   //    Newly allocated profile object
 161   //
 162   TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
 163                                    Form("ESD/MC signal for FMD%d%c", d, r),
 164                                    axis.GetNbins(),
 165                                    axis.GetXmin(),
 166                                    axis.GetXmax(),
 167                                    (r == 'I' || r == 'i') ? 20 : 40,
 168                                    0, 2*TMath::Pi());
 169   ret->GetXaxis()->SetTitle("#eta");
 170   ret->GetYaxis()->SetTitle("#varphi [radians]");
 171   ret->GetZaxis()->SetTitle("#LT #deltaN_{ch,incl}=N_{ch,incl,ana}-N_{ch,incl,mc}#GT");
 172   ret->SetDirectory(0);
 173   return ret;
 174 }
 175 //____________________________________________________________________
 176 TH2D*
 177 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
 178 {
 179   //
 180   // MAke comparison profiles
 181   //
 182   // Parameters:
 183   //    d     Detector
 184   //    r     Ring
 185   //
 186   // Return:
 187   //    Newly allocated profile object
 188   //
 189   TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
 190                        Form("ESD-MC correlation for FMD%d%c", d, r),
 191                        41, -.5, 40.5, 410, -.5, 40.5);
 192   ret->GetXaxis()->SetTitle("m_{incl} (MC)");
 193   ret->GetYaxis()->SetTitle("m_{incl} from #Delta/#Delta_{mp} (ESD)");
 194   ret->GetZaxis()->SetTitle("Correlation");
 195   ret->SetDirectory(0);
 196   return ret;
 197 }
 198 //____________________________________________________________________
 199 TH1D*
 200 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r, Int_t max) const
 201 {
 202   //
 203   // MAke comparison profiles
 204   //
 205   // Parameters:
 206   //    d     Detector
 207   //    r     Ring
 208   //
 209   // Return:
 210   //    Newly allocated profile object
 211   //
 212   TH1D* ret = new TH1D(Form("FMD%d%c_diff_mc_esd", d, r),
 213                        Form("ESD-MC difference for FMD%d%c", d, r),
 214                        5*max,-max, max);
 215   ret->GetXaxis()->SetTitle("MC - ESD");
 216   ret->SetFillColor(AliForwardUtil::RingColor(d,r));
 217   ret->SetLineColor(AliForwardUtil::RingColor(d,r));
 218   ret->SetFillStyle(3001);
 219   ret->SetDirectory(0);
 220   return ret;
 221 }
 222 //____________________________________________________________________
 223 void
 224 AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
 225 {
 226   //
 227   // Fill comparison profiles
 228   //
 229   // Parameters:
 230   //    d    Detector
 231   //    r    Ring
 232   //    esd  ESD histogram
 233   //    mc   MC histogram
 234   //
 235   if (!esd || !mc) return;
 236   TProfile2D* p = 0;
 237   TH2D*       h = 0;
 238   TH1D*       x = 0;
 239   switch (d) {
 240   case 1:
 241     p = fFMD1i;
 242     h = fFMD1iC;
 243     x = fFMD1iD;
 244     break;
 245   case 2:
 246     p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o);
 247     h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC);
 248     x = (r == 'I' || r == 'i' ? fFMD2iD : fFMD2oD);
 249     break;
 250   case 3:
 251     p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o);
 252     h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC);
 253     x = (r == 'I' || r == 'i' ? fFMD3iD : fFMD3oD);
 254     break;
 255   }
 256   if (!p) return;
 257
 258   for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
 259     Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
 260     for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
 261       Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
 262       Double_t mEsd = esd->GetBinContent(iEta,iPhi);
 263       Double_t mMc  = mc->GetBinContent(iEta,iPhi);
 264
 265       if (mMc <= 0 && mEsd <= 0) continue;
 266       x->Fill(mMc-mEsd);
 267       p->Fill(eta, phi, mEsd - mMc);
 268       h->Fill(mMc,mEsd);
 269     }
 270   }
 271 }
 272
 273 //____________________________________________________________________
 274 Bool_t
 275 AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
 276                                           AliForwardUtil::Histos& mc)
 277 {
 278   //
 279   // Compare the result of analysing the ESD for
 280   // the inclusive charged particle density to analysing
 281   // MC truth
 282   //
 283   // Parameters:
 284   //    esd
 285   //    mc
 286   //
 287   // Return:
 288   //   true
 289   //
 290   Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
 291   Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
 292   Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
 293   Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
 294   Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
 295
 296   return kTRUE;
 297 }
 298
 299 //____________________________________________________________________
 300 void
 301 AliFMDMCDensityCalculator::CreateOutputObjects(TList* dir)
 302 {
 303   //
 304   // Output diagnostic histograms to directory
 305   //
 306   // Parameters:
 307   //    dir List to write in
 308   //
 309   AliFMDDensityCalculator::CreateOutputObjects(dir);
 310   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 311
 312   fComps = new TList;
 313   fComps->SetOwner();
 314   fComps->SetName("esd_mc_comparison");
 315   d->Add(fComps);
 316
 317 }
 318 //____________________________________________________________________
 319 //
 320 // EOF
 321 //
 322
 323
 324