PWGLF/FORWARD/analysis2/AliForwarddNdetaTask.cxx

   1 //====================================================================
   2 #include "AliForwarddNdetaTask.h"
   3 #include <TMath.h>
   4 #include <TH2D.h>
   5 #include <TH1D.h>
   6 #include <THStack.h>
   7 #include <TList.h>
   8 #include <TFile.h>
   9 #include <AliAnalysisManager.h>
  10 #include <AliAODEvent.h>
  11 #include <AliAODHandler.h>
  12 #include <AliAODInputHandler.h>
  13 #include "AliForwardUtil.h"
  14 #include "AliAODForwardMult.h"
  15
  16 //____________________________________________________________________
  17 AliForwarddNdetaTask::AliForwarddNdetaTask()
  18   : AliBasedNdetaTask()
  19 {
  20   //
  21   // Constructor
  22   //
  23   DGUARD(fDebug, 3, "Default CTOR of AliForwarddNdetaTask");
  24 }
  25
  26 //____________________________________________________________________
  27 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
  28   : AliBasedNdetaTask("Forward")
  29 {
  30   //
  31   // Constructor
  32   //
  33   // Paramters
  34   //   name    Name of task
  35   SetTitle("FMD");
  36   DGUARD(fDebug, 3, "Named CTOR of AliForwarddNdetaTask");
  37 }
  38
  39 //____________________________________________________________________
  40 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
  41   : AliBasedNdetaTask(o)
  42 {
  43   //
  44   // Copy constructor
  45   //
  46   DGUARD(fDebug, 3, "Copy CTOR of AliForwarddNdetaTask");
  47 }
  48
  49 //____________________________________________________________________
  50 AliBasedNdetaTask::CentralityBin*
  51 AliForwarddNdetaTask::MakeCentralityBin(const char* name, Short_t l, Short_t h)
  52   const
  53 {
  54   //
  55   // Make a new centrality bin
  56   //
  57   // Parameters:
  58   //    name   Histogram names
  59   //    l      Lower cut
  60   //    h      Upper cut
  61   //
  62   // Return:
  63   //    Newly allocated object (of our type)
  64   //
  65   DGUARD(fDebug, 3,"Make a centrality bin for AliForwarddNdetaTask: %s [%d,%d]",
  66          name, l, h);
  67   return new AliForwarddNdetaTask::CentralityBin(name, l, h);
  68 }
  69
  70
  71 //____________________________________________________________________
  72 TH2D*
  73 AliForwarddNdetaTask::GetHistogram(const AliAODEvent* aod, Bool_t mc)
  74 {
  75   //
  76   // Retrieve the histogram
  77   //
  78   // Parameters:
  79   //    aod AOD event
  80   //    mc  Whether to get the MC histogram or not
  81   //
  82   // Return:
  83   //    Retrieved histogram or null
  84   //
  85   TObject* obj = 0;
  86   if (mc) obj = aod->FindListObject("ForwardMC");
  87   else    obj = aod->FindListObject("Forward");
  88
  89   // We should have a forward object at least
  90   if (!obj) {
  91     if (!mc) AliWarning("No Forward object found AOD");
  92     return 0;
  93   }
  94   AliAODForwardMult* forward = static_cast<AliAODForwardMult*>(obj);
  95   return &(forward->GetHistogram());
  96 }
  97 //____________________________________________________________________
  98 void
  99 AliForwarddNdetaTask::CheckEventData(Double_t vtx,
 100                                      TH2*     data,
 101                                      TH2*     dataMC)
 102 {
 103   // Check if this is satellite
 104   // if (!fSatelliteVertices) return;
 105   Double_t aVtx = TMath::Abs(vtx);
 106   if (aVtx < 37.5 || aVtx > 400) return;
 107
 108   TH2* hists[] = { data, dataMC };
 109
 110   // In satellite vertices FMD2i is cut away manually at this point
 111   // for certain vertices. It could be done in the ESDs, but as of
 112   // this writing not for specific vertices.
 113   //
 114   // cholm comment: It would be difficult to setup the filter in the
 115   // reconstruction pass, but it could perhaps be done in the AOD
 116   // filtering.
 117   //
 118   // This is what was done for
 119   // the Pb-Pb paper (arXiv:1304.0347).
 120   for (Int_t iX = 0; iX<=data->GetNbinsX(); iX++) {
 121     // Do all checks up front - as soon as we can - branching is
 122     // expensive!
 123     Double_t x    = data->GetXaxis()->GetBinCenter(iX);
 124     Bool_t   zero = false;
 125     if (((vtx >  60 && vtx <  90) && x < 3) ||
 126         ((vtx > 330 && vtx < 350) && x > -2.5) ||
 127         ((vtx < 100 || vtx > 305) && TMath::Abs(x) < 4.5) ||
 128         (vtx < 50                 && TMath::Abs(x) < 4.75))
 129       zero = true;
 130     if (!zero) continue;
 131
 132     for (Int_t iH = 0; iH < 2; iH++) {
 133       if (!hists[iH]) continue;
 134       // if (iX > hists[iH]->GetNbinsX()+1) continue;
 135       // Also zero coverage and phi acceptance for this
 136       for (Int_t iY = 0; iY<=hists[iH]->GetNbinsY()+1; iY++) {
 137         hists[iH]->SetBinContent(iX, iY, 0);
 138         hists[iH]->SetBinError(iX, iY, 0);
 139       }
 140     }
 141   }
 142
 143   if (fCorrEmpty) {
 144     // Now, since we have some dead areas in FMD2i (sectors 16 and
 145     // 17), we need to remove the corresponding bins from the
 146     // histogram. However, it is not obvious which bins (in eta) to
 147     // remove, so remove everything starting from the most negative to
 148     // the middle of the histogram.
 149     //
 150     // This hack was first introduced by HHD, but was done at the end of
 151     // the event processing (CentralityBin::MakeResults).  That is,
 152     // however, not very practical, as we'd like to normalize to the phi
 153     // acceptance rather than the eta coverage and then correct for
 154     // empty bins. Since the only way to really update the phi
 155     // acceptance stored in the overflow bin is on the event level, we
 156     // should really do it here.
 157     const Int_t phiBin1 = 17; // Sector 16
 158     const Int_t phiBin2 = 18; // Sector 17
 159     for (Int_t iH = 0; iH < 2; iH++) {
 160       if (!hists[iH]) continue;
 161
 162       Int_t midX = hists[iH]->GetNbinsX() / 2;
 163       // Int_t nY   = hists[iH]->GetNbinsY();
 164       for (Int_t i = 1; i <= midX; i++) {
 165         hists[iH]->SetBinContent(i, phiBin1, 0);
 166         hists[iH]->SetBinContent(i, phiBin2, 0);
 167         hists[iH]->SetBinError(i, phiBin1, 0);
 168         hists[iH]->SetBinError(i, phiBin2, 0);
 169
 170         // Here, we should also modify the overflow bin to reflect the
 171         // new phi acceptance.  First get the old phi acceptance -
 172         // then multiply this on the number of bins. This gives us -
 173         // roughly - the number of sectors we had.  Then take out two
 174         // from that number, and then calculate the new phi
 175         // Acceptance. Note, if the sectors where already taken out in
 176         // the AOD production, we _will_ end up with a wrong number,
 177         // so we should _not_ do that in the AOD production.  This is
 178         // tricky and may not work at all.  For now, we should rely on
 179         // the old way of correcting to the eta coverage and
 180         // correcting for empty bins.
 181       }
 182     }
 183   }
 184 }
 185 //========================================================================
 186 void
 187 AliForwarddNdetaTask::CentralityBin::End(TList*      sums,
 188                                          TList*      results,
 189                                          UShort_t    scheme,
 190                                          const TH2F* shapeCorr,
 191                                          Double_t    trigEff,
 192                                          Double_t    trigEff0,
 193                                          Bool_t      symmetrice,
 194                                          Int_t       rebin,
 195                                          Bool_t      rootProj,
 196                                          Bool_t      corrEmpty,
 197                                          Bool_t      cutEdges,
 198                                          Int_t       triggerMask,
 199                                          Int_t       marker,
 200                                          Int_t       color,
 201                                          TList*      mclist,
 202                                          TList*      truthlist )
 203 {
 204   DGUARD(fDebug, 1, "In End of %s with corrEmpty=%d, cutEdges=%d, rootProj=%d",
 205          GetName(), corrEmpty, cutEdges, rootProj);
 206   AliBasedNdetaTask::CentralityBin::End(sums, results, scheme,
 207                                         shapeCorr, trigEff, trigEff0,
 208                                         symmetrice, rebin,
 209                                         rootProj, corrEmpty, cutEdges,
 210                                         triggerMask, marker, color, mclist,
 211                                         truthlist);
 212
 213   if (!IsAllBin()) return;
 214   TFile* file = TFile::Open("forward.root", "READ");
 215   if (!file) return;
 216
 217   TList* forward = static_cast<TList*>(file->Get("Forward"));
 218   if (!forward) {
 219     AliError("List Forward not found in forward.root");
 220     return;
 221   }
 222   TList* rings = static_cast<TList*>(forward->FindObject("ringResults"));
 223   if (!rings) {
 224     AliError("List ringResults not found in forward.root");
 225     return;
 226   }
 227   THStack* res = static_cast<THStack*>(rings->FindObject("all"));
 228   if (!res) {
 229     AliError(Form("Stack all not found in %s", rings->GetName()));
 230     return;
 231   }
 232   if (!fTriggers) {
 233     AliError("Triggers histogram not set");
 234     return;
 235   }
 236
 237   Double_t ntotal   = 0;
 238   Double_t epsilonT = trigEff;
 239 #if 0
 240   // TEMPORARY FIX
 241   if (triggerMask == AliAODForwardMult::kNSD) {
 242     // This is a local change
 243     epsilonT = 0.92;
 244     AliWarning(Form("Using hard-coded NSD trigger efficiency of %f",epsilonT));
 245   }
 246 #endif
 247   AliInfo("Adding per-ring histograms to output");
 248   TString text;
 249   Double_t scaler = Normalization(*fTriggers, scheme, epsilonT, ntotal, &text);
 250   TIter next(res->GetHists());
 251   TH1*  hist = 0;
 252   while ((hist = static_cast<TH1*>(next()))) hist->Scale(scaler);
 253   res->SetName("dndetaRings");
 254   fOutput->Add(res);
 255   fOutput->Add(new TNamed("normCalc", text.Data()));
 256 }
 257
 258 //________________________________________________________________________
 259 //
 260 // EOF
 261 //