PWG2/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->Clear();
  37   if (fFMD1i)  delete fFMD1i;
  38   if (fFMD2i)  delete fFMD2i;
  39   if (fFMD2o)  delete fFMD2o;
  40   if (fFMD3i)  delete fFMD3i;
  41   if (fFMD3o)  delete fFMD3o;
  42   if (fFMD1iC) delete fFMD1iC;
  43   if (fFMD2iC) delete fFMD2iC;
  44   if (fFMD2oC) delete fFMD2oC;
  45   if (fFMD3iC) delete fFMD3iC;
  46   if (fFMD3oC) delete fFMD3oC;
  47 }
  48
  49 //____________________________________________________________________
  50 AliFMDMCDensityCalculator&
  51 AliFMDMCDensityCalculator::operator=(const AliFMDMCDensityCalculator& o)
  52 {
  53   //
  54   // Assignement operator
  55   //
  56   // Parameters:
  57   //    o Object to assign from
  58   //
  59   // Return:
  60   //    Reference to this object
  61   //
  62   AliFMDDensityCalculator::operator=(o);
  63   return *this;
  64 }
  65
  66
  67 //____________________________________________________________________
  68 void
  69 AliFMDMCDensityCalculator::Init(const TAxis& eAxis)
  70 {
  71   //
  72   // Initialize this object
  73   //
  74   // Parameters:
  75   //    etaAxis Eta axis to use
  76   //
  77   AliFMDDensityCalculator::Init(eAxis);
  78   fFMD1i  = Make(1,'I',eAxis);
  79   fFMD2i  = Make(2,'I',eAxis);
  80   fFMD2o  = Make(2,'O',eAxis);
  81   fFMD3i  = Make(3,'I',eAxis);
  82   fFMD3o  = Make(3,'O',eAxis);
  83   fFMD1iC = Make(1,'I');
  84   fFMD2iC = Make(2,'I');
  85   fFMD2oC = Make(2,'O');
  86   fFMD3iC = Make(3,'I');
  87   fFMD3oC = Make(3,'O');
  88
  89   fComps->Add(fFMD1i);
  90   fComps->Add(fFMD2i);
  91   fComps->Add(fFMD2o);
  92   fComps->Add(fFMD3i);
  93   fComps->Add(fFMD3o);
  94   fComps->Add(fFMD1iC);
  95   fComps->Add(fFMD2iC);
  96   fComps->Add(fFMD2oC);
  97   fComps->Add(fFMD3iC);
  98   fComps->Add(fFMD3oC);
  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   //    event  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, s, t, eta, 0, 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 [degrees]");
 171   ret->GetZaxis()->SetTitle("#LT incusive density ESD/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                        200, 0, 20, 200, 0, 20);
 192   ret->GetXaxis()->SetTitle("m_{incl} (MC)");
 193   ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
 194   ret->GetZaxis()->SetTitle("Correlation");
 195   ret->SetDirectory(0);
 196   return ret;
 197 }
 198 //____________________________________________________________________
 199 void
 200 AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
 201 {
 202   //
 203   // Fill comparison profiles
 204   //
 205   // Parameters:
 206   //    d    Detector
 207   //    r    Ring
 208   //    esd  ESD histogram
 209   //    mc   MC histogram
 210   //
 211   if (!esd || !mc) return;
 212   TProfile2D* p = 0;
 213   TH2D*       h = 0;
 214   switch (d) {
 215   case 1:
 216     p = fFMD1i;
 217     h = fFMD1iC;
 218     break;
 219   case 2:
 220     p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o);
 221     h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC);
 222     break;
 223   case 3:
 224     p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o);
 225     h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC);
 226     break;
 227   }
 228   if (!p) return;
 229
 230   for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
 231     Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
 232     for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
 233       Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
 234       Double_t mEsd = esd->GetBinContent(iEta,iPhi);
 235       Double_t mMc  = mc->GetBinContent(iEta,iPhi);
 236
 237       p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
 238       h->Fill(mMc,mEsd);
 239     }
 240   }
 241 }
 242
 243 //____________________________________________________________________
 244 Bool_t
 245 AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
 246                                           AliForwardUtil::Histos& mc)
 247 {
 248   //
 249   // Compare the result of analysing the ESD for
 250   // the inclusive charged particle density to analysing
 251   // MC truth
 252   //
 253   // Parameters:
 254   //    esd
 255   //    mc
 256   //
 257   // Return:
 258   //   true
 259   //
 260   Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
 261   Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
 262   Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
 263   Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
 264   Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
 265
 266   return kTRUE;
 267 }
 268
 269 //____________________________________________________________________
 270 void
 271 AliFMDMCDensityCalculator::DefineOutput(TList* dir)
 272 {
 273   //
 274   // Output diagnostic histograms to directory
 275   //
 276   // Parameters:
 277   //    dir List to write in
 278   //
 279   AliFMDDensityCalculator::DefineOutput(dir);
 280   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 281
 282   fComps = new TList;
 283   fComps->SetName("esd_mc_comparison");
 284   d->Add(fComps);
 285
 286 }
 287 //____________________________________________________________________
 288 //
 289 // EOF
 290 //
 291
 292
 293