PWG2/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 <AliAnalysisManager.h>
   9 #include <AliAODEvent.h>
  10 #include <AliAODHandler.h>
  11 #include <AliAODInputHandler.h>
  12 #include "AliForwardUtil.h"
  13
  14 //____________________________________________________________________
  15 AliForwarddNdetaTask::AliForwarddNdetaTask()
  16   : AliAnalysisTaskSE(),
  17     fSumForward(0),     //  Sum of histograms
  18     fSumForwardMC(0),   //  Sum of MC histograms (if any)
  19     fSumPrimary(0),     //  Sum of primary histograms
  20     fSumCentral(0),     //  Sum of central histograms
  21     fCentral(0),        // Cache of central histogram
  22     fPrimary(0),        // Cache of primary histogram
  23     fSums(0),           // Container of sums
  24     fOutput(0),         // Container of outputs
  25     fTriggers(0),       // Histogram of triggers
  26     fVtxMin(0),         // Minimum v_z
  27     fVtxMax(0),         // Maximum v_z
  28     fTriggerMask(0),    // Trigger mask
  29     fRebin(0),          // Rebinning factor
  30     fCutEdges(false),
  31     fSNNString(0),
  32     fSysString(0)
  33 {}
  34
  35 //____________________________________________________________________
  36 AliForwarddNdetaTask::AliForwarddNdetaTask(const char* /* name */)
  37   : AliAnalysisTaskSE("Forward"),
  38     fSumForward(0),     //  Sum of histograms
  39     fSumForwardMC(0),   //  Sum of MC histograms (if any)
  40     fSumPrimary(0),     //  Sum of primary histograms
  41     fSumCentral(0),     //  Sum of central histograms
  42     fCentral(0),        // Cache of central histogram
  43     fPrimary(0),        // Cache of primary histogram
  44     fSums(0),           // Container of sums
  45     fOutput(0),         // Container of outputs
  46     fTriggers(0),       // Histogram of triggers
  47     fVtxMin(-10),       // Minimum v_z
  48     fVtxMax(10),        // Maximum v_z
  49     fTriggerMask(AliAODForwardMult::kInel),
  50     fRebin(5),          // Rebinning factor
  51     fCutEdges(false),
  52     fSNNString(0),
  53     fSysString(0)
  54 {
  55   // Output slot #1 writes into a TH1 container
  56   DefineOutput(1, TList::Class());
  57   DefineOutput(2, TList::Class());
  58 }
  59
  60 //____________________________________________________________________
  61 AliForwarddNdetaTask::AliForwarddNdetaTask(const AliForwarddNdetaTask& o)
  62   : AliAnalysisTaskSE(o),
  63     fSumForward(o.fSumForward),         // TH2D* -  Sum of histograms
  64     fSumForwardMC(o.fSumForwardMC),             // TH2D* -  Sum of MC histograms (if any)
  65     fSumPrimary(o.fSumPrimary),         // TH2D* -  Sum of primary histograms
  66     fSumCentral(o.fSumCentral),         // TH2D* -  Sum of central histograms
  67     fCentral(o.fCentral),       //! Cache of central histogram
  68     fPrimary(o.fPrimary),       //! Cache of primary histogram
  69     fSums(o.fSums),             // TList* - Container of sums
  70     fOutput(o.fOutput),         // TList* - Container of outputs
  71     fTriggers(o.fTriggers),             // TH1D* - Histogram of triggers
  72     fVtxMin(o.fVtxMin),         // Double_t - Minimum v_z
  73     fVtxMax(o.fVtxMax),         // Double_t - Maximum v_z
  74     fTriggerMask(o.fTriggerMask),               // Int_t - Trigger mask
  75     fRebin(o.fRebin),           // Int_t - Rebinning factor
  76     fCutEdges(o.fCutEdges),             // Bool_t - Whether to cut edges when rebinning
  77     fSNNString(o.fSNNString),           // TNamed* -
  78     fSysString(o.fSysString)            // TNamed* -
  79 {}
  80
  81 //____________________________________________________________________
  82 AliForwarddNdetaTask::~AliForwarddNdetaTask()
  83 {
  84   if (fSums) {
  85     fSums->Delete();
  86     delete fSums;
  87     fSums = 0;
  88   }
  89   if (fOutputs) {
  90     fOutputs->Delete();
  91     delete fOutputs;
  92     fOutputs = 0;
  93   }
  94   if (fCentral)    delete fCentral;
  95   if (fPrimary)    delete fPrimary;
  96 }
  97
  98 //________________________________________________________________________
  99 void
 100 AliForwarddNdetaTask::SetTriggerMask(const char* mask)
 101 {
 102   UShort_t    trgMask = 0;
 103   TString     trgs(mask);
 104   trgs.ToUpper();
 105   TObjString* trg;
 106   TIter       next(trgs.Tokenize(" ,|"));
 107   while ((trg = static_cast<TObjString*>(next()))) {
 108     TString s(trg->GetString());
 109     if      (s.IsNull()) continue;
 110     if      (s.CompareTo("INEL")  == 0) trgMask = AliAODForwardMult::kInel;
 111     else if (s.CompareTo("INEL>0")== 0) trgMask = AliAODForwardMult::kInelGt0;
 112     else if (s.CompareTo("NSD")   == 0) trgMask = AliAODForwardMult::kNSD;
 113     else
 114       Warning("SetTriggerMask", "Unknown trigger %s", s.Data());
 115   }
 116   if (trgMask == 0) trgMask = 1;
 117   SetTriggerMask(trgMask);
 118 }
 119
 120 //________________________________________________________________________
 121 void
 122 AliForwarddNdetaTask::UserCreateOutputObjects()
 123 {
 124   // Create histograms
 125   // Called once (on the worker node)
 126
 127   fOutput = new TList;
 128   fOutput->SetName(Form("%s_result", GetName()));
 129   fOutput->SetOwner();
 130
 131   fSums = new TList;
 132   fSums->SetName(Form("%s_sums", GetName()));
 133   fSums->SetOwner();
 134
 135
 136   fTriggers = new TH1D("triggers", "Number of triggers",
 137                        kAccepted, 1, kAccepted);
 138   fTriggers->SetYTitle("# of events");
 139   fTriggers->GetXaxis()->SetBinLabel(kAll,         "All events");
 140   fTriggers->GetXaxis()->SetBinLabel(kB,           "w/B trigger");
 141   fTriggers->GetXaxis()->SetBinLabel(kA,           "w/A trigger");
 142   fTriggers->GetXaxis()->SetBinLabel(kC,           "w/C trigger");
 143   fTriggers->GetXaxis()->SetBinLabel(kE,           "w/E trigger");
 144   fTriggers->GetXaxis()->SetBinLabel(kMB,          "w/Collision trigger");
 145   fTriggers->GetXaxis()->SetBinLabel(kWithVertex,  "w/Vertex");
 146   fTriggers->GetXaxis()->SetBinLabel(kWithTrigger, "w/Selected trigger");
 147   fTriggers->GetXaxis()->SetBinLabel(kAccepted,    "Accepted by cut");
 148   fTriggers->GetXaxis()->SetNdivisions(kAccepted, false);
 149   fTriggers->SetFillColor(kRed+1);
 150   fTriggers->SetFillStyle(3001);
 151   fTriggers->SetStats(0);
 152   fSums->Add(fTriggers);
 153
 154   fSNNString = new TNamed("sNN", "");
 155   fSysString = new TNamed("sys", "");
 156   fSums->Add(fSNNString);
 157   fSums->Add(fSysString);
 158
 159   // Check that we have an AOD input handler
 160   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
 161   AliAODInputHandler* ah =
 162     dynamic_cast<AliAODInputHandler*>(am->GetInputEventHandler());
 163   if (!ah) AliFatal("No AOD input handler set in analysis manager");
 164
 165   // Post data for ALL output slots >0 here, to get at least an empty histogram
 166   PostData(1, fSums);
 167 }
 168
 169 //____________________________________________________________________
 170 TH2D*
 171 AliForwarddNdetaTask::CloneHist(TH2D* in, const char* name)
 172 {
 173   if (!in) return 0;
 174   TH2D* ret = static_cast<TH2D*>(in->Clone(name));
 175   ret->SetDirectory(0);
 176   ret->Sumw2();
 177   fSums->Add(ret);
 178
 179   return ret;
 180 }
 181
 182 //____________________________________________________________________
 183 void
 184 AliForwarddNdetaTask::UserExec(Option_t *)
 185 {
 186   // Main loop
 187   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(InputEvent());
 188   if (!aod) {
 189     AliError("Cannot get the AOD event");
 190     return;
 191   }
 192
 193   // Get objects from the event structure
 194   TObject* oForward   = aod->FindListObject("Forward");
 195   TObject* oForwardMC = aod->FindListObject("ForwardMC");
 196   TObject* oPrimary   = aod->FindListObject("primary");
 197   TObject* oCentral   = aod->FindListObject("Central");
 198
 199   // We should have a forward object at least
 200   if (!oForward) {
 201     AliWarning("No Forward object found AOD");
 202     return;
 203   }
 204
 205   // Cast to good types
 206   AliAODForwardMult* forward   = static_cast<AliAODForwardMult*>(oForward);
 207   AliAODForwardMult* forwardMC = static_cast<AliAODForwardMult*>(oForwardMC);
 208   TH2D*              primary   = static_cast<TH2D*>(oPrimary);
 209   TH2D*              central   = static_cast<TH2D*>(oCentral);
 210
 211   static bool first = true;
 212   if (first) {
 213     UShort_t sNN = forward->GetSNN();
 214     UShort_t sys = forward->GetSystem();
 215     fSNNString->SetTitle(AliForwardUtil::CenterOfMassEnergyString(sNN));
 216     fSNNString->SetUniqueID(sNN);
 217     fSysString->SetTitle(AliForwardUtil::CollisionSystemString(sys));
 218     fSysString->SetUniqueID(sys);
 219     first = false;
 220   }
 221
 222   // Create our sum histograms
 223   if (!fSumForward) fSumForward = CloneHist(&forward->GetHistogram(),"forward");
 224   if (!fSumForwardMC && forwardMC)
 225     fSumForwardMC = CloneHist(&forwardMC->GetHistogram(),"forwardMC");
 226   if (!fSumPrimary && primary) fSumPrimary = CloneHist(primary, "primary");
 227   if (!fSumCentral && central) fSumCentral = CloneHist(central, "central");
 228
 229   // Add contribtion from MC
 230   if (primary) fSumPrimary->Add(primary);
 231
 232   // Count event
 233   fTriggers->AddBinContent(kAll);
 234   if (forward->IsTriggerBits(AliAODForwardMult::kB))
 235     fTriggers->AddBinContent(kB);
 236   if (forward->IsTriggerBits(AliAODForwardMult::kA))
 237     fTriggers->AddBinContent(kA);
 238   if (forward->IsTriggerBits(AliAODForwardMult::kC))
 239     fTriggers->AddBinContent(kC);
 240   if (forward->IsTriggerBits(AliAODForwardMult::kE))
 241     fTriggers->AddBinContent(kE);
 242   if (forward->IsTriggerBits(AliAODForwardMult::kInel))
 243     fTriggers->AddBinContent(kMB);
 244
 245   // Check if we have an event of interest.
 246   if (!forward->IsTriggerBits(fTriggerMask)) return;
 247   fTriggers->AddBinContent(kWithTrigger);
 248
 249   // Check that we have a valid vertex
 250   if (!forward->HasIpZ()) return;
 251   fTriggers->AddBinContent(kWithVertex);
 252
 253   // Check that vertex is within cuts
 254   if (!forward->InRange(fVtxMin, fVtxMax)) return;
 255   fTriggers->AddBinContent(kAccepted);
 256
 257   // Add contribution
 258   fSumForward->Add(&(forward->GetHistogram()));
 259   if (fSumForwardMC) fSumForwardMC->Add(&(forwardMC->GetHistogram()));
 260   if (fSumPrimary)   fSumPrimary->Add(primary);
 261   if (fSumCentral)   fSumCentral->Add(central);
 262
 263   PostData(1, fSums);
 264 }
 265
 266 //________________________________________________________________________
 267 void
 268 AliForwarddNdetaTask::SetHistogramAttributes(TH1D* h, Int_t colour, Int_t marker,
 269                                           const char* title, const char* ytitle)
 270 {
 271   h->SetTitle(title);
 272   h->SetMarkerColor(colour);
 273   h->SetMarkerStyle(marker);
 274   h->SetMarkerSize(1);
 275   h->SetFillStyle(0);
 276   h->SetYTitle(ytitle);
 277   h->GetXaxis()->SetTitleFont(132);
 278   h->GetXaxis()->SetLabelFont(132);
 279   h->GetXaxis()->SetNdivisions(10);
 280   h->GetYaxis()->SetTitleFont(132);
 281   h->GetYaxis()->SetLabelFont(132);
 282   h->GetYaxis()->SetNdivisions(10);
 283   h->GetYaxis()->SetDecimals();
 284   h->SetStats(0);
 285 }
 286
 287 //________________________________________________________________________
 288 void
 289 AliForwarddNdetaTask::Terminate(Option_t *)
 290 {
 291   // Draw result to screen, or perform fitting, normalizations Called
 292   // once at the end of the query
 293
 294   fSums = dynamic_cast<TList*> (GetOutputData(1));
 295   if(!fOutput) {
 296     AliError("Could not retrieve TList fSums");
 297     return;
 298   }
 299
 300   if (!fOutput) {
 301     fOutput = new TList;
 302     fOutput->SetName(Form("%s_result", GetName()));
 303     fOutput->SetOwner();
 304   }
 305
 306   fSumForward     = static_cast<TH2D*>(fSums->FindObject("forward"));
 307   fSumForwardMC   = static_cast<TH2D*>(fSums->FindObject("forwardMC"));
 308   fSumPrimary     = static_cast<TH2D*>(fSums->FindObject("primary"));
 309   fSumCentral     = static_cast<TH2D*>(fSums->FindObject("central"));
 310   fTriggers       = static_cast<TH1D*>(fSums->FindObject("triggers"));
 311   fSNNString      = static_cast<TNamed*>(fSums->FindObject("sNN"));
 312   fSysString      = static_cast<TNamed*>(fSums->FindObject("sys"));
 313
 314   if (!fTriggers) {
 315     AliError("Couldn't find histogram 'triggers' in list");
 316     return;
 317   }
 318   if (!fSumForward) {
 319     AliError("Couldn't find histogram 'forward' in list");
 320     return;
 321   }
 322
 323   Int_t nAll        = Int_t(fTriggers->GetBinContent(kAll));
 324   Int_t nB          = Int_t(fTriggers->GetBinContent(kB));
 325   Int_t nA          = Int_t(fTriggers->GetBinContent(kA));
 326   Int_t nC          = Int_t(fTriggers->GetBinContent(kC));
 327   Int_t nE          = Int_t(fTriggers->GetBinContent(kE));
 328   Int_t nMB         = Int_t(fTriggers->GetBinContent(kMB));
 329   Int_t nTriggered  = Int_t(fTriggers->GetBinContent(kWithTrigger));
 330   Int_t nWithVertex = Int_t(fTriggers->GetBinContent(kWithVertex));
 331   Int_t nAccepted   = Int_t(fTriggers->GetBinContent(kAccepted));
 332   Int_t nGood       = nB - nA - nC + 2 * nE;
 333   Double_t vtxEff   = Double_t(nMB) / nTriggered * Double_t(nAccepted) / nGood;
 334   AliInfo(Form("Total of %9d events\n"
 335                "                   of these %9d are minimum bias\n"
 336                "                   of these %9d has a %s trigger\n"
 337                "                   of these %9d has a vertex\n"
 338                "                   of these %9d were in [%+4.1f,%+4.1f]cm\n"
 339                "                   Triggers by type:\n"
 340                "                     B   = %9d\n"
 341                "                     A|C = %9d (%9d+%-9d)\n"
 342                "                     E   = %9d\n"
 343                "                   Implies %9d good triggers\n"
 344                "                   Vertex efficiency: %f",
 345                nAll, nMB, nTriggered,
 346                AliAODForwardMult::GetTriggerString(fTriggerMask),
 347                nWithVertex, nAccepted,
 348                fVtxMin, fVtxMax,
 349                nB, nA+nC, nA, nC, nE, nGood, vtxEff));
 350
 351   // Get acceptance normalisation from underflow bins
 352   TH1D* normForward   = fSumForward->ProjectionX("normForward", 0, 1, "");
 353   // Project onto eta axis - _ignoring_underflow_bins_!
 354   TH1D* dndetaForward = fSumForward->ProjectionX("dndetaForward", 1, -1, "e");
 355   // Normalize to the acceptance
 356   dndetaForward->Divide(normForward);
 357   // Scale by the vertex efficiency
 358   dndetaForward->Scale(vtxEff, "width");
 359
 360   SetHistogramAttributes(dndetaForward, kRed+1, 20, "ALICE Forward");
 361   SetHistogramAttributes(normForward, kRed+1, 20, "ALICE Forward", "Normalisation");
 362
 363   fOutput->Add(fTriggers->Clone());
 364   fOutput->Add(fSNNString->Clone());
 365   fOutput->Add(fSysString->Clone());
 366   fOutput->Add(dndetaForward);
 367   fOutput->Add(normForward);
 368   fOutput->Add(Rebin(dndetaForward));
 369
 370   if (fSumForwardMC) {
 371     // Get acceptance normalisation from underflow bins
 372     TH1D* normForwardMC   = fSumForwardMC->ProjectionX("normForwardMC", 0, 1, "");
 373     // Project onto eta axis - _ignoring_underflow_bins_!
 374     TH1D* dndetaForwardMC = fSumForwardMC->ProjectionX("dndetaForwardMC", 1, -1, "e");
 375     // Normalize to the acceptance
 376     dndetaForwardMC->Divide(normForwardMC);
 377     // Scale by the vertex efficiency
 378     dndetaForwardMC->Scale(vtxEff, "width");
 379
 380     SetHistogramAttributes(dndetaForwardMC, kRed+3, 21, "ALICE Forward (MC)");
 381     SetHistogramAttributes(normForwardMC, kRed+3, 21, "ALICE Forward (MC)",
 382                            "Normalisation");
 383
 384     fOutput->Add(dndetaForwardMC);
 385     fOutput->Add(normForwardMC);
 386     fOutput->Add(Rebin(dndetaForwardMC));
 387   }
 388
 389   if (fSumPrimary) {
 390     TH1D* dndetaTruth = fSumPrimary->ProjectionX("dndetaTruth",-1,-1,"e");
 391     dndetaTruth->Scale(1./nAll, "width");
 392
 393     SetHistogramAttributes(dndetaTruth, kGray+3, 22, "Monte-Carlo truth");
 394
 395     fOutput->Add(dndetaTruth);
 396     fOutput->Add(Rebin(dndetaTruth));
 397   }
 398   if (fSumCentral) {
 399     TH1D* dndetaCentral = fSumCentral->ProjectionX("dndetaCentral",-1,-1,"e");
 400     dndetaCentral->Scale(1./nGood, "width");
 401
 402     SetHistogramAttributes(dndetaCentral, kGray+3, 22, "ALICE Central - track(let)s");
 403
 404     dndetaCentral->GetXaxis()->SetRangeUser(-1,1);
 405     fOutput->Add(dndetaCentral);
 406     fOutput->Add(Rebin(dndetaCentral));
 407   }
 408
 409   TNamed* trigString =
 410     new TNamed("trigString", AliAODForwardMult::GetTriggerString(fTriggerMask));
 411   trigString->SetUniqueID(fTriggerMask);
 412   fOutput->Add(trigString);
 413
 414   TAxis* vtxAxis = new TAxis(1,fVtxMin,fVtxMax);
 415   vtxAxis->SetName("vtxAxis");
 416   vtxAxis->SetTitle(Form("v_{z}#in[%+5.1f,%+5.1f]cm", fVtxMin,fVtxMax));
 417   fOutput->Add(vtxAxis);
 418
 419   // If only there was a away to get sqrt{s_NN} and beam type
 420
 421
 422   PostData(2, fOutput);
 423 }
 424
 425 //________________________________________________________________________
 426 TH1D*
 427 AliForwarddNdetaTask::Rebin(const TH1D* h) const
 428 {
 429   if (fRebin <= 1) return 0;
 430
 431   Int_t nBins = h->GetNbinsX();
 432   if(nBins % fRebin != 0) {
 433     Warning("Rebin", "Rebin factor %d is not a devisor of current number "
 434             "of bins %d in the histogram %s", fRebin, nBins, h->GetName());
 435     return 0;
 436   }
 437
 438   // Make a copy
 439   TH1D* tmp = static_cast<TH1D*>(h->Clone(Form("%s_rebin%02d",
 440                                                h->GetName(), fRebin)));
 441   tmp->Rebin(fRebin);
 442   tmp->SetDirectory(0);
 443
 444   // The new number of bins
 445   Int_t nBinsNew = nBins / fRebin;
 446   for(Int_t i = 1;i<= nBinsNew; i++) {
 447     Double_t content = 0;
 448     Double_t sumw    = 0;
 449     Double_t wsum    = 0;
 450     Int_t    nbins   = 0;
 451     for(Int_t j = 1; j<=fRebin;j++) {
 452       Int_t    bin = (i-1)*fRebin + j;
 453       Double_t c   =  h->GetBinContent(bin);
 454
 455       if (c <= 0) continue;
 456
 457       if (fCutEdges) {
 458         if (h->GetBinContent(bin+1)<=0 ||
 459             h->GetBinContent(bin-1)) {
 460           Warning("Rebin", "removing bin %d=%f of %s (%d=%f,%d=%f)",
 461                   bin, c, h->GetName(),
 462                   bin+1, h->GetBinContent(bin+1),
 463                   bin-1, h->GetBinContent(bin-1));
 464           continue;
 465         }
 466       }
 467       Double_t e =  h->GetBinError(bin);
 468       Double_t w =  1 / (e*e); // 1/c/c
 469       content    += c;
 470       sumw       += w;
 471       wsum       += w * c;
 472       nbins++;
 473     }
 474
 475     if(content > 0 && nbins > 0) {
 476       tmp->SetBinContent(i, wsum / sumw);
 477       tmp->SetBinError(i,1./TMath::Sqrt(sumw));
 478     }
 479   }
 480
 481   return tmp;
 482 }