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