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   fFMD1i  = Make(1,'I',eAxis);
  78   fFMD2i  = Make(2,'I',eAxis);
  79   fFMD2o  = Make(2,'O',eAxis);
  80   fFMD3i  = Make(3,'I',eAxis);
  81   fFMD3o  = Make(3,'O',eAxis);
  82   fFMD1iC = Make(1,'I');
  83   fFMD2iC = Make(2,'I');
  84   fFMD2oC = Make(2,'O');
  85   fFMD3iC = Make(3,'I');
  86   fFMD3oC = Make(3,'O');
  87
  88   fComps->Add(fFMD1i);
  89   fComps->Add(fFMD2i);
  90   fComps->Add(fFMD2o);
  91   fComps->Add(fFMD3i);
  92   fComps->Add(fFMD3o);
  93   fComps->Add(fFMD1iC);
  94   fComps->Add(fFMD2iC);
  95   fComps->Add(fFMD2oC);
  96   fComps->Add(fFMD3iC);
  97   fComps->Add(fFMD3oC);
  98 }
  99
 100
 101 //____________________________________________________________________
 102 Bool_t
 103 AliFMDMCDensityCalculator::CalculateMC(const AliESDFMD&        fmd,
 104                                        AliForwardUtil::Histos& hists)
 105 {
 106   //
 107   // Calculate the charged particle density from the MC track references.
 108   //
 109   // Parameters:
 110   //    event  MC event
 111   //    hists  Histograms to fill
 112   //    vz     Interaction z coordinate @f$ v_z@f$
 113   //    vtxBin bin corresponding to @f$ v_z@f$
 114   //
 115   // Return:
 116   //    true on success
 117   //
 118   for (UShort_t d=1; d<=3; d++) {
 119     UShort_t nr = (d == 1 ? 1 : 2);
 120     for (UShort_t q=0; q<nr; q++) {
 121       Char_t      r = (q == 0 ? 'I' : 'O');
 122       UShort_t    ns= (q == 0 ?  20 :  40);
 123       UShort_t    nt= (q == 0 ? 512 : 256);
 124       TH2D*       h = hists.Get(d,r);
 125
 126       for (UShort_t s=0; s<ns; s++) {
 127         for (UShort_t t=0; t<nt; t++) {
 128           Float_t mult = fmd.Multiplicity(d,r,s,t);
 129
 130           if (mult == 0 || mult > 20) continue;
 131
 132           Float_t phi = fmd.Phi(d,r,s,t) / 180 * TMath::Pi();
 133           Float_t eta = fmd.Eta(d,r,s,t);
 134
 135           Float_t corr = Correction(d, r, s, t, eta, 0, false);
 136           Float_t sig  = (corr <= 0 ? 0 : mult / corr);
 137           h->Fill(eta,phi,sig);
 138         }
 139       }
 140     }
 141   }
 142   return kTRUE;
 143 }
 144
 145 //____________________________________________________________________
 146 TProfile2D*
 147 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r,
 148                                 const TAxis& axis) const
 149 {
 150   //
 151   // MAke comparison profiles
 152   //
 153   // Parameters:
 154   //    d     Detector
 155   //    r     Ring
 156   //    axis  Eta axis
 157   //
 158   // Return:
 159   //    Newly allocated profile object
 160   //
 161   TProfile2D* ret = new TProfile2D(Form("FMD%d%c_esd_vs_mc", d, r),
 162                                    Form("ESD/MC signal for FMD%d%c", d, r),
 163                                    axis.GetNbins(),
 164                                    axis.GetXmin(),
 165                                    axis.GetXmax(),
 166                                    (r == 'I' || r == 'i') ? 20 : 40,
 167                                    0, 2*TMath::Pi());
 168   ret->GetXaxis()->SetTitle("#eta");
 169   ret->GetYaxis()->SetTitle("#varphi [degrees]");
 170   ret->GetZaxis()->SetTitle("#LT incusive density ESD/MC#GT");
 171   ret->SetDirectory(0);
 172   return ret;
 173 }
 174 //____________________________________________________________________
 175 TH2D*
 176 AliFMDMCDensityCalculator::Make(UShort_t d, Char_t r) const
 177 {
 178   //
 179   // MAke comparison profiles
 180   //
 181   // Parameters:
 182   //    d     Detector
 183   //    r     Ring
 184   //
 185   // Return:
 186   //    Newly allocated profile object
 187   //
 188   TH2D* ret = new TH2D(Form("FMD%d%c_corr_mc_esd", d, r),
 189                        Form("ESD-MC correlation for FMD%d%c", d, r),
 190                        200, 0, 20, 200, 0, 20);
 191   ret->GetXaxis()->SetTitle("m_{incl} (MC)");
 192   ret->GetYaxis()->SetTitle("#Delta/#Delta_{mp} (ESD)");
 193   ret->GetZaxis()->SetTitle("Correlation");
 194   ret->SetDirectory(0);
 195   return ret;
 196 }
 197 //____________________________________________________________________
 198 void
 199 AliFMDMCDensityCalculator::Fill(UShort_t d, Char_t r, TH2* esd, TH2* mc)
 200 {
 201   //
 202   // Fill comparison profiles
 203   //
 204   // Parameters:
 205   //    d    Detector
 206   //    r    Ring
 207   //    esd  ESD histogram
 208   //    mc   MC histogram
 209   //
 210   if (!esd || !mc) return;
 211   TProfile2D* p = 0;
 212   TH2D*       h = 0;
 213   switch (d) {
 214   case 1:
 215     p = fFMD1i;
 216     h = fFMD1iC;
 217     break;
 218   case 2:
 219     p = (r == 'I' || r == 'i' ? fFMD2i  : fFMD2o);
 220     h = (r == 'I' || r == 'i' ? fFMD2iC : fFMD2oC);
 221     break;
 222   case 3:
 223     p = (r == 'I' || r == 'i' ? fFMD3i  : fFMD3o);
 224     h = (r == 'I' || r == 'i' ? fFMD3iC : fFMD3oC);
 225     break;
 226   }
 227   if (!p) return;
 228
 229   for (Int_t iEta = 1; iEta <= esd->GetNbinsX(); iEta++) {
 230     Double_t eta = esd->GetXaxis()->GetBinCenter(iEta);
 231     for (Int_t iPhi = 1; iPhi <= esd->GetNbinsY(); iPhi++) {
 232       Double_t phi  = esd->GetYaxis()->GetBinCenter(iPhi);
 233       Double_t mEsd = esd->GetBinContent(iEta,iPhi);
 234       Double_t mMc  = mc->GetBinContent(iEta,iPhi);
 235
 236       p->Fill(eta, phi, (mMc > 0 ? mEsd / mMc : 0));
 237       h->Fill(mMc,mEsd);
 238     }
 239   }
 240 }
 241
 242 //____________________________________________________________________
 243 Bool_t
 244 AliFMDMCDensityCalculator::CompareResults(AliForwardUtil::Histos& esd,
 245                                           AliForwardUtil::Histos& mc)
 246 {
 247   //
 248   // Compare the result of analysing the ESD for
 249   // the inclusive charged particle density to analysing
 250   // MC truth
 251   //
 252   // Parameters:
 253   //    esd
 254   //    mc
 255   //
 256   // Return:
 257   //   true
 258   //
 259   Fill(1, 'I', esd.Get(1,'I'), mc.Get(1,'I'));
 260   Fill(2, 'I', esd.Get(2,'I'), mc.Get(2,'I'));
 261   Fill(2, 'O', esd.Get(2,'O'), mc.Get(2,'O'));
 262   Fill(3, 'I', esd.Get(3,'I'), mc.Get(3,'I'));
 263   Fill(3, 'O', esd.Get(3,'O'), mc.Get(3,'O'));
 264
 265   return kTRUE;
 266 }
 267
 268 //____________________________________________________________________
 269 void
 270 AliFMDMCDensityCalculator::DefineOutput(TList* dir)
 271 {
 272   //
 273   // Output diagnostic histograms to directory
 274   //
 275   // Parameters:
 276   //    dir List to write in
 277   //
 278   AliFMDDensityCalculator::DefineOutput(dir);
 279   TList* d = static_cast<TList*>(dir->FindObject(GetName()));
 280
 281   fComps = new TList;
 282   fComps->SetName("esd_mc_comparison");
 283   d->Add(fComps);
 284
 285 }
 286 //____________________________________________________________________
 287 //
 288 // EOF
 289 //
 290
 291
 292